Integration across different omics datasets

1 Introduction

[Needs more text]

In this document we present a work-flow for integration across different omics datasets.

[Note] This is not the final version of the document.

2 Packages and files

A package with regularized CCA and multiomics DIABLO method is mixOmics. Package igraph is needed for network analysis.

library(mixOmics)
library(igraph)

Package for complex heatmaps.

# BiocManager::install('ComplexHeatmap')
library(ComplexHeatmap)
# https://www.rdocumentation.org/packages/pheatmap/versions/1.0.12/topics/pheatmap
library(pheatmap)

2.1 Functions

Some original and adapted functions can be found in the file that is silently processed here.

out <- ""
out <- paste(out, knit_child("005-Functions.Rmd", quiet = TRUE))

2.2 File names

Usually we use a file management system under the pISA-tree framework. For simplicity, all files ( scripts, data ,… ) are in one directory.

Sample description file (aka phenodata)

pfn <- "../input/phenodata_subset_2023-03-08.txt"
files <- c("../input/data_Proteomics.txt", "../input/data_metabolomics.txt", 
    "../input/data_hormonomics.txt", "../input/data_qPCR.txt", "../input/data_Phenomics.txt")

For future use and labeling, we need text names of dataset objects.

datanames <- c("proteomics", "metabolomics", "hormonomics", "qPCR", "phenomics")
datanames

## [1] "proteomics"   "metabolomics" "hormonomics"  "qPCR"        
## [5] "phenomics"

Define groups to consider. In a small example, we will pick two groups, based on treatment:

useTreatment <- c("H", "HD")

Just to check: Phenodata file name

pfn

## [1] "../input/phenodata_subset_2023-03-08.txt"

Data file names

files

## [1] "../input/data_Proteomics.txt"   "../input/data_metabolomics.txt"
## [3] "../input/data_hormonomics.txt"  "../input/data_qPCR.txt"        
## [5] "../input/data_Phenomics.txt"

It is advisable to first read the phenodata, followed by actual data input. This enables smart selection of samples, based on the sample selection column with the assay name.

2.3 Phenodata

pfn

## [1] "../input/phenodata_subset_2023-03-08.txt"

# 
phdata <- read.table(pfn, header = TRUE, sep = "\t", stringsAsFactors = FALSE, 
    row.names = 1)
dim(phdata)

## [1] 56 30

names(phdata)

##  [1] "Treatment"                        
##  [2] "Harvest"                          
##  [3] "SamplingDay"                      
##  [4] "DaysOfStressH"                    
##  [5] "DaysOfStressD"                    
##  [6] "DaysOfStressW"                    
##  [7] "DaysRecovery"                     
##  [8] "Replicate"                        
##  [9] "Sample.type"                      
## [10] "Date"                             
## [11] "Heat.Drought.Water.Recovery.Days" 
## [12] "TreatmentxDatexPlant"             
## [13] "TreatmentxSamplingDay"            
## [14] "TreatmentxSamplingDayxReplicateNo"
## [15] "Comment"                          
## [16] "Num_Tubers"                       
## [17] "Total_Tubers_Weight"              
## [18] "FW_SH_Total"                      
## [19] "DW_SH_Total"                      
## [20] "Fv.Fm"                            
## [21] "qL"                               
## [22] "deltaT"                           
## [23] "TOP.AREA"                         
## [24] "COMPACTNESS"                      
## [25] "WATER.CONSUMPTION"                
## [26] "Proteomics"                       
## [27] "Transcriptomics"                  
## [28] "Metabolomics"                     
## [29] "Hormonomics"                      
## [30] "X_A_300_OverView.R"

pdata <- phdata

Filter out groups

pdata$Treatment

##  [1] "C"  "C"  "C"  "C"  "C"  "C"  "C"  "C"  "C"  "C"  "C"  "C"  "C" 
## [14] "C"  "C"  "C"  "D"  "D"  "D"  "D"  "D"  "D"  "D"  "D"  "H"  "H" 
## [27] "H"  "H"  "H"  "H"  "H"  "H"  "H"  "H"  "H"  "H"  "H"  "H"  "H" 
## [40] "H"  "HD" "HD" "HD" "HD" "HD" "HD" "HD" "HD" "W"  "W"  "W"  "W" 
## [53] "W"  "W"  "W"  "W"

pdata <- pdata[pdata$Treatment %in% useTreatment, ]
pdata$Treatment

##  [1] "H"  "H"  "H"  "H"  "H"  "H"  "H"  "H"  "H"  "H"  "H"  "H"  "H" 
## [14] "H"  "H"  "H"  "HD" "HD" "HD" "HD" "HD" "HD" "HD" "HD"

dim(pdata)

## [1] 24 30

Overview of selected samples:

addmargins(table(pdata$Treatment, pdata$SamplingDay))

##      
##        1  7  8 14 Sum
##   H    4  4  4  4  16
##   HD   0  0  4  4   8
##   Sum  4  4  8  8  24

.treat <- unique(pdata$Treatment)
.days <- unique(pdata$SamplingDay)
.entry <- 0.5

summary(pdata)

##   Treatment            Harvest       SamplingDay     DaysOfStressH   
##  Length:24          Min.   :1.000   Min.   : 1.000   Min.   : 1.000  
##  Class :character   1st Qu.:2.000   1st Qu.: 7.000   1st Qu.: 7.000  
##  Mode  :character   Median :3.000   Median : 8.000   Median : 8.000  
##                     Mean   :2.833   Mean   : 8.667   Mean   : 8.667  
##                     3rd Qu.:4.000   3rd Qu.:14.000   3rd Qu.:14.000  
##                     Max.   :4.000   Max.   :14.000   Max.   :14.000  
##  DaysOfStressD   DaysOfStressW  DaysRecovery   Replicate   
##  Min.   :0.000   Min.   :0     Min.   :0     Min.   :1.00  
##  1st Qu.:0.000   1st Qu.:0     1st Qu.:0     1st Qu.:1.75  
##  Median :0.000   Median :0     Median :0     Median :2.50  
##  Mean   :1.333   Mean   :0     Mean   :0     Mean   :2.50  
##  3rd Qu.:1.000   3rd Qu.:0     3rd Qu.:0     3rd Qu.:3.25  
##  Max.   :7.000   Max.   :0     Max.   :0     Max.   :4.00  
##  Sample.type            Date          
##  Length:24          Length:24         
##  Class :character   Class :character  
##  Mode  :character   Mode  :character  
##                                       
##                                       
##                                       
##  Heat.Drought.Water.Recovery.Days TreatmentxDatexPlant
##  Length:24                        Length:24           
##  Class :character                 Class :character    
##  Mode  :character                 Mode  :character    
##                                                       
##                                                       
##                                                       
##  TreatmentxSamplingDay TreatmentxSamplingDayxReplicateNo
##  Length:24             Length:24                        
##  Class :character      Class :character                 
##  Mode  :character      Mode  :character                 
##                                                         
##                                                         
##                                                         
##    Comment          Num_Tubers     Total_Tubers_Weight  FW_SH_Total
##  Length:24          Mode:logical   Mode:logical        Min.   :1   
##  Class :character   NA's:24        NA's:24             1st Qu.:1   
##  Mode  :character                                      Median :1   
##                                                        Mean   :1   
##                                                        3rd Qu.:1   
##                                                        Max.   :1   
##   DW_SH_Total     Fv.Fm         qL        deltaT     TOP.AREA
##  Min.   :1    Min.   :1   Min.   :1   Min.   :1   Min.   :1  
##  1st Qu.:1    1st Qu.:1   1st Qu.:1   1st Qu.:1   1st Qu.:1  
##  Median :1    Median :1   Median :1   Median :1   Median :1  
##  Mean   :1    Mean   :1   Mean   :1   Mean   :1   Mean   :1  
##  3rd Qu.:1    3rd Qu.:1   3rd Qu.:1   3rd Qu.:1   3rd Qu.:1  
##  Max.   :1    Max.   :1   Max.   :1   Max.   :1   Max.   :1  
##   COMPACTNESS WATER.CONSUMPTION   Proteomics Transcriptomics
##  Min.   :1    Min.   :1         Min.   :1    Min.   :1      
##  1st Qu.:1    1st Qu.:1         1st Qu.:1    1st Qu.:1      
##  Median :1    Median :1         Median :1    Median :1      
##  Mean   :1    Mean   :1         Mean   :1    Mean   :1      
##  3rd Qu.:1    3rd Qu.:1         3rd Qu.:1    3rd Qu.:1      
##  Max.   :1    Max.   :1         Max.   :1    Max.   :1      
##   Metabolomics  Hormonomics X_A_300_OverView.R
##  Min.   :1     Min.   :1    Min.   : 1.000    
##  1st Qu.:1     1st Qu.:1    1st Qu.: 7.000    
##  Median :1     Median :1    Median : 8.000    
##  Mean   :1     Mean   :1    Mean   : 8.667    
##  3rd Qu.:1     3rd Qu.:1    3rd Qu.:14.000    
##  Max.   :1     Max.   :1    Max.   :14.000

apply(pdata, 2, function(x) summary(as.factor(x)))

## $Treatment
##  H HD 
## 16  8 
## 
## $Harvest
## 1 2 3 4 
## 4 4 8 8 
## 
## $SamplingDay
##  1  7  8 14 
##  4  4  8  8 
## 
## $DaysOfStressH
##  1  7  8 14 
##  4  4  8  8 
## 
## $DaysOfStressD
##  0  1  7 
## 16  4  4 
## 
## $DaysOfStressW
##  0 
## 24 
## 
## $DaysRecovery
##  0 
## 24 
## 
## $Replicate
## 1 2 3 4 
## 6 6 6 6 
## 
## $Sample.type
##  S 
## 24 
## 
## $Date
## 2020-11-04 2020-11-10 2020-11-11 2020-11-17 
##          4          4          8          8 
## 
## $Heat.Drought.Water.Recovery.Days
##  1_0_0_0 14_0_0_0 14_7_0_0  7_0_0_0  8_0_0_0  8_1_0_0 
##        4        4        4        4        4        4 
## 
## $TreatmentxDatexPlant
##  H_2020-11-04_1  H_2020-11-04_2  H_2020-11-04_3  H_2020-11-04_4 
##               1               1               1               1 
##  H_2020-11-10_1  H_2020-11-10_2  H_2020-11-10_3  H_2020-11-10_4 
##               1               1               1               1 
##  H_2020-11-11_1  H_2020-11-11_2  H_2020-11-11_3  H_2020-11-11_4 
##               1               1               1               1 
##  H_2020-11-17_1  H_2020-11-17_2  H_2020-11-17_3  H_2020-11-17_4 
##               1               1               1               1 
## HD_2020-11-11_1 HD_2020-11-11_2 HD_2020-11-11_3 HD_2020-11-11_4 
##               1               1               1               1 
## HD_2020-11-17_1 HD_2020-11-17_2 HD_2020-11-17_3 HD_2020-11-17_4 
##               1               1               1               1 
## 
## $TreatmentxSamplingDay
##   H_1  H_14   H_7   H_8 HD_14  HD_8 
##     4     4     4     4     4     4 
## 
## $TreatmentxSamplingDayxReplicateNo
##   H_1_1   H_1_2   H_1_3   H_1_4  H_14_1  H_14_2  H_14_3  H_14_4 
##       1       1       1       1       1       1       1       1 
##   H_7_1   H_7_2   H_7_3   H_7_4   H_8_1   H_8_2   H_8_3   H_8_4 
##       1       1       1       1       1       1       1       1 
## HD_14_1 HD_14_2 HD_14_3 HD_14_4  HD_8_1  HD_8_2  HD_8_3  HD_8_4 
##       1       1       1       1       1       1       1       1 
## 
## $Comment
## only for cor 
##           24 
## 
## $Num_Tubers
## NA's 
##   24 
## 
## $Total_Tubers_Weight
## NA's 
##   24 
## 
## $FW_SH_Total
##  1 
## 24 
## 
## $DW_SH_Total
##  1 
## 24 
## 
## $Fv.Fm
##  1 
## 24 
## 
## $qL
##  1 
## 24 
## 
## $deltaT
##  1 
## 24 
## 
## $TOP.AREA
##  1 
## 24 
## 
## $COMPACTNESS
##  1 
## 24 
## 
## $WATER.CONSUMPTION
##  1 
## 24 
## 
## $Proteomics
##  1 
## 24 
## 
## $Transcriptomics
##  1 
## 24 
## 
## $Metabolomics
##  1 
## 24 
## 
## $Hormonomics
##  1 
## 24 
## 
## $X_A_300_OverView.R
##  1  7  8 14 
##  4  4  8  8

2.4 Process data files

For this project we aim to integrate several multi-omics datasets. We have data on hormonomics, metabolomics, and qPCR:

Read all datafiles and assign them to named objects

for (i in 1:length(files)) {
    print(datanames[i])
    assign(datanames[i], read.table(files[i], header = TRUE, sep = "\t"))
    print(head(get(datanames[i])))
}

## [1] "proteomics"
##   SampleID VdnDe2_86784 PBdnRY1_7025 VdnPW5_1541 CLCdnPW10_11437
## 1   C_S1_1  0.000413884  0.000469004 0.002308758              NA
## 2   C_S1_2  0.000230027  0.000469191 0.002352451     0.000300619
## 3   C_S1_3  0.000449712  0.000524164 0.001314037              NA
## 4   C_S1_4  0.000255355  0.000520852 0.001661846     0.000222480
## 5   C_S7_1  0.000178758           NA 0.000498582     0.000622978
## 6   C_S7_2  0.000226049  0.000307384 0.000910694     0.000787786
##   CLCdnDe2_96512 VdnPW1_81435 CLCdnDe6_4953 CLCdnDe13_5738
## 1             NA     2.82e-05   0.000386853             NA
## 2             NA     2.35e-05   0.000322506             NA
## 3             NA     2.62e-05   0.000520422             NA
## 4             NA           NA   0.000397796             NA
## 5             NA           NA   0.000222778             NA
## 6             NA           NA   0.000563427             NA
##   CLCdnPW12_1472 CLCdnPW49_8560 VdnPW4_18801 SdnRY1_10712
## 1    0.000767766    0.000039000  0.000203446  0.000272489
## 2    0.000474884             NA  0.000197874  0.000265026
## 3    0.000530523    0.000072600  0.000189478  0.000338374
## 4    0.000584473    0.000036100  0.000188281  0.000336236
## 5    0.000288815    0.000101007  0.000439347  0.000470757
## 6    0.000432853             NA  0.000407422  0.000347256
##   CLCdnDe2_311 CLCdnDe13_63179 CLCdnDe10_24187 TDe_comp110944_c2_seq4
## 1  0.000257351     0.000135326     0.000119484            0.000050600
## 2  0.000178787     0.000037600     0.000033200                     NA
## 3  0.000159788     0.000042000     0.000037100            0.000094200
## 4  0.000238167     0.000041700     0.000073700            0.000046800
## 5  0.000222302     0.000116896              NA            0.000131029
## 6  0.000187408     0.000246367              NA            0.000165693
##   PBdnRY1_11085 CLCdnPW48_15164 CLCdnPW31_956 VdnDe1_162036
## 1      4.44e-05     0.000263690   0.000334185   0.000146521
## 2      1.85e-05     0.000302265   0.000167159   0.000122150
## 3      4.14e-05     0.000276283   0.000145245   0.000091000
## 4            NA     0.000274538   0.000247419   0.000180799
## 5            NA     0.000427083   0.000288672   0.000253132
## 6            NA     0.000216027   0.000243360   0.000160049
##   PBdnRY1_2239 CLCdnPW22_15452 CLCdnPW50_2470 VdnPW3_37011 SdnPW1_46
## 1  0.000165383     0.000161016             NA  0.000123402  1.49e-05
## 2  0.000068900     0.000089500             NA  0.000085700        NA
## 3  0.000057800     0.000100000             NA  0.000057500        NA
## 4  0.000057400     0.000149014             NA  0.000114203  1.38e-05
## 5  0.000107145              NA             NA  0.000106596        NA
## 6  0.000045200     0.000058600             NA  0.000112329  3.26e-05
##   CLCdnDe7_34340 CLCdnDe11_2286 VdnDe2_53978 CLCdnDe7_6144
## 1             NA             NA  0.000323764   0.000176082
## 2             NA             NA  0.000337389   0.000366984
## 3             NA             NA  0.000452302   0.000245989
## 4             NA             NA  0.000524353   0.000244435
## 5             NA    0.000369451           NA   0.000228152
## 6             NA    0.000311459  0.001149384   0.000192340
##   CLCdnRY10_6741 PBdnRY1_5389 VdnDe2_29290 VdnDe4_115100
## 1    0.000141861  0.000708110  0.000141561   0.000711383
## 2    0.000059100  0.000641660  0.000094400   0.000773551
## 3    0.000066100  0.000602145  0.000158210   0.000547316
## 4    0.000065600  0.000598341  0.000183412   0.000572483
## 5             NA  0.000877616  0.000073400   0.001041977
## 6    0.000077500  0.000807118  0.000092800   0.001047347
##   Sotub11g025210.1.1 VdnDe1_39992 VdnPW4_15664
## 1        0.000228539  0.000389774           NA
## 2        0.000285788  0.000541569           NA
## 3        0.000319272  0.000423515           NA
## 4        0.000211504  0.000360720           NA
## 5        0.000197414  0.000505036           NA
## 6        0.000332853  0.000425761           NA
## [1] "metabolomics"
##   SampleID Glukose Fructose Sucrose Starch    Asp    Glu   Asn   Ser
## 1   C_S1_1    2.13     2.70   3.449  22.06 1039.5 2514.0 177.8 597.9
## 2   C_S1_2    2.20     2.90   3.382  12.74  844.3 1966.4 167.8 498.1
## 3   C_S1_3    0.82     1.59   2.452   9.98  887.2 2067.8 167.0 441.1
## 4   C_S1_4    2.55     3.01   4.057  15.13  987.6 2347.9 172.3 537.6
## 5   C_S7_1    4.77     3.97   3.533  16.25  793.2 2109.0 276.9 368.5
## 6   C_S7_2    6.30     7.16   6.280   9.57 1391.6 3344.4 498.8 704.2
##     Gln   Gly  His  Arg   Thr   Ala   Pro  Tyr  Val Met  Ile  Lys
## 1 497.8 137.8 20.7 26.9 225.4 702.5  48.6 25.3 54.1 7.3 48.8 26.0
## 2 408.8 104.6 13.3 29.1 208.5 495.5  57.3 22.1 51.3 6.9 47.0 21.8
## 3 400.1  92.1 17.1 29.5 216.1 514.6  53.5 24.4 52.9 8.6 54.3 26.5
## 4 465.7 117.2 16.8 24.9 252.0 653.1  58.6 22.3 58.4 7.9 52.9 26.7
## 5 265.7  68.7 13.7 38.3 188.6 296.4  85.8 31.2 73.1 4.7 60.4 28.4
## 6 680.9 101.2 19.4 63.9 302.4 664.6 135.7 40.6 96.3 8.6 80.6 33.1
##    Leu   Phe
## 1 18.8 119.0
## 2 15.7  88.0
## 3 19.5  86.8
## 4 20.1 111.0
## 5 27.4  98.3
## 6 27.1 171.2
## [1] "hormonomics"
##   SampleID      IAA    oxIAA  IAA.Asp      ABA        PA       DPA
## 1   C_S1_1 37.19565 61.49305 2.224072 36.78497  91.99991  45.63410
## 2   C_S1_2 45.86649 67.84222 1.990994 33.11733  92.61196  62.10651
## 3   C_S1_3 47.60516 52.49759 1.496438 41.66091  93.83119  55.06159
## 4   C_S1_4 37.55562 48.69611 2.240201 39.11995  91.35024  59.82576
## 5   C_S7_1 49.46230 76.79460 1.927162 80.00281 231.68757 178.08837
## 6   C_S7_2 78.91379 93.14059 1.856353 49.76028 166.08573 149.86757
##          SA       JA    JA.Ile X9.10.dhJA X12.OH.JA   cisOPDA
## 1  505.2129 2.687303 0.5525212   5.449105  16.13451  353.8608
## 2  519.3793 2.802456 0.5661384   3.700174  23.95458  402.7049
## 3  275.2972 4.761132 0.4273195   2.879761  27.93628  644.9688
## 4  628.0462 4.622662 0.2028412   5.165337  25.72584  750.3532
## 5 1314.9488 6.101443 0.6226329   5.010268 244.12270 1294.9082
## 6  731.2122 3.092070 0.3449961   8.373161 203.74167  873.4217
## [1] "qPCR"
##   SampleID    RbohA    SnRK2     ACO2    HSP70     PR1b    RD29B
## 1   C_S1_1 1.346213 1.497998 0.150379 1.013451 0.324007 0.016995
## 2   C_S1_2 1.496014 1.627412 0.196017 1.067230 0.154344 0.046504
## 3   C_S1_3 1.262812 1.630533 0.440406 0.883357 0.269488 0.016995
## 4   C_S1_4 1.428252 1.530426 0.176572 1.042068 0.255647 0.082387
## 5   C_S7_1 1.404150 1.122600 0.537281 1.028170 0.226714 1.751881
## 6   C_S7_2 0.848707 0.725752 0.150857 0.640827 0.192341 0.113202
##    X13.LOX     P5CS     ERF1     CAT1       CO    SWEET     SP6A
## 1 0.699873 3.647878 0.560424 0.639920 2.421594 0.816298 0.122374
## 2 0.660413 3.264675 0.637808 0.635107 5.563090 1.873813 0.238811
## 3 0.765764 2.151917 0.585558 0.687279 2.740914 0.930861 0.330151
## 4 0.734167 3.155558 0.663618 0.704880 2.469193 0.933900 0.122374
## 5 1.136201 0.862709 1.634810 0.811006 1.709219 1.494808 3.386457
## 6 0.540666 0.908966 1.239857 0.373886 1.580358 0.853875 0.810631
##     M0ZJG3
## 1 1.211115
## 2 1.376397
## 3 1.007479
## 4 0.903495
## 5 2.357865
## 6 1.387190
## [1] "phenomics"
##   SampleID FW_SH_Total DW_SH_Total  TOP.AREA COMPACTNESS     Fv.Fm
## 1   C_S1_1      62.511       2.939  90507.25   0.7167274 0.7065431
## 2   C_S1_2      83.262       4.388 101872.30   0.7121689 0.7133399
## 3   C_S1_3      65.553       3.199  92754.84   0.7108707 0.7137180
## 4   C_S1_4      73.533       3.808 108909.81   0.6882440 0.7064169
## 5   C_S7_1     112.044       7.998 156269.47   0.7066436 0.6953883
## 6   C_S7_2     101.276       6.851 154364.77   0.6910989 0.6940050
##          qL    deltaT WATER.CONSUMPTION
## 1 0.8462101 1.3552322                77
## 2 0.8415086 1.2458286                68
## 3 0.7952288 0.8838253                78
## 4 0.8433725 1.0340176                94
## 5 0.8848586 1.2201157               132
## 6 0.9179832 1.0551987               133

Check if all 5 objects are created

datanames

## [1] "proteomics"   "metabolomics" "hormonomics"  "qPCR"        
## [5] "phenomics"

datanames %in% ls()

## [1] TRUE TRUE TRUE TRUE TRUE

2.5 Combine datasets

Datasets for DIABLO need to be collected in a list, with rows corresponding to the same samples. The order of samples from shrinked phenodata will be enforced.

The first component of the list will be a grouping variable, indicating the conditions. We will create reasonable names for groups.

sday <- paste0(0, pdata$SamplingDay)
len <- nchar(sday)
sday <- substr(sday, len - 1, len)
trt <- pdata$Treatment
what <- paste(trt, sday, sep = "")
what

##  [1] "H01"  "H01"  "H01"  "H01"  "H07"  "H07"  "H07"  "H07"  "H08" 
## [10] "H08"  "H08"  "H08"  "H14"  "H14"  "H14"  "H14"  "HD08" "HD08"
## [19] "HD08" "HD08" "HD14" "HD14" "HD14" "HD14"

X <- list(status = what)
names(X[[1]]) <- rownames(pdata)
X

## $status
##   H_S1_1   H_S1_2   H_S1_3   H_S1_4   H_S7_1   H_S7_2   H_S7_3 
##    "H01"    "H01"    "H01"    "H01"    "H07"    "H07"    "H07" 
##   H_S7_4   H_S8_1   H_S8_2   H_S8_3   H_S8_4  H_S14_1  H_S14_2 
##    "H07"    "H08"    "H08"    "H08"    "H08"    "H14"    "H14" 
##  H_S14_3  H_S14_4  HD_S8_1  HD_S8_2  HD_S8_3  HD_S8_4 HD_S14_1 
##    "H14"    "H14"   "HD08"   "HD08"   "HD08"   "HD08"   "HD14" 
## HD_S14_2 HD_S14_3 HD_S14_4 
##   "HD14"   "HD14"   "HD14"

print(addmargins(table(pdata$SamplingDay, what)), zero.print = ".")

##      what
##       H01 H07 H08 H14 HD08 HD14 Sum
##   1     4   .   .   .    .    .   4
##   7     .   4   .   .    .    .   4
##   8     .   .   4   .    4    .   8
##   14    .   .   .   4    .    4   8
##   Sum   4   4   4   4    4    4  24

print(addmargins(table(pdata$Treatment, what)), zero.print = ".")

##      what
##       H01 H07 H08 H14 HD08 HD14 Sum
##   H     4   4   4   4    .    .  16
##   HD    .   .   .   .    4    4   8
##   Sum   4   4   4   4    4    4  24

Put datasets into the list X and ensure that they all have same order of samples as in phenodata.

datanames

## [1] "proteomics"   "metabolomics" "hormonomics"  "qPCR"        
## [5] "phenomics"

i <- 1
for (i in 1:length(datanames)) {
    x <- get(datanames[i])
    rownames(x) <- x[, 1]
    x <- x[, -1]
    X[[i + 1]] <- x[rownames(pdata), ]
    names(X)[i + 1] <- datanames[i]
}
str(X)

## List of 6
##  $ status      : Named chr [1:24] "H01" "H01" "H01" "H01" ...
##   ..- attr(*, "names")= chr [1:24] "H_S1_1" "H_S1_2" "H_S1_3" "H_S1_4" ...
##  $ proteomics  :'data.frame':    24 obs. of  36 variables:
##   ..$ VdnDe2_86784          : num [1:24] 0.000281 0.000522 0.00038 0.000467 0.000374 ...
##   ..$ PBdnRY1_7025          : num [1:24] 0.000478 0.000443 0.000603 0.000529 0.000407 ...
##   ..$ VdnPW5_1541           : num [1:24] 0.00131 0.00218 0.00212 0.00152 0.00104 ...
##   ..$ CLCdnPW10_11437       : num [1:24] 0.000612 0.000455 0.000441 0.000543 0.000782 ...
##   ..$ CLCdnDe2_96512        : num [1:24] NA NA NA NA NA ...
##   ..$ VdnPW1_81435          : num [1:24] 2.87e-05 2.66e-05 2.58e-05 3.18e-05 NA NA NA NA NA NA ...
##   ..$ CLCdnDe6_4953         : num [1:24] 0.000569 0.000447 0.000395 0.000485 0.000512 ...
##   ..$ CLCdnDe13_5738        : num [1:24] NA NA NA NA 0.000149 ...
##   ..$ CLCdnPW12_1472        : num [1:24] 0.000542 0.000503 0.000466 0.000545 0.000483 ...
##   ..$ CLCdnPW49_8560        : num [1:24] 3.97e-05 3.69e-05 NA NA 8.45e-05 ...
##   ..$ VdnPW4_18801          : num [1:24] 0.000276 0.000192 0.000218 0.000268 0.000404 ...
##   ..$ SdnRY1_10712          : num [1:24] 0.00037 0.000386 0.000375 0.000359 0.000295 ...
##   ..$ CLCdnDe2_311          : num [1:24] 0.000262 0.000122 0.000157 0.000194 0.000186 ...
##   ..$ CLCdnDe13_63179       : num [1:24] 9.19e-05 8.53e-05 4.14e-05 2.55e-04 1.96e-04 ...
##   ..$ CLCdnDe10_24187       : num [1:24] NA NA NA NA 0.00013 ...
##   ..$ TDe_comp110944_c2_seq4: num [1:24] 5.15e-05 4.78e-05 NA 5.71e-05 1.64e-04 ...
##   ..$ PBdnRY1_11085         : num [1:24] NA NA NA NA NA NA NA NA 6.20e-05 7.43e-05 ...
##   ..$ CLCdnPW48_15164       : num [1:24] 0.000201 0.00028 0.000333 0.000186 0.000357 ...
##   ..$ CLCdnPW31_956         : num [1:24] 0.000136 0.00019 0.000184 0.000251 0.000217 ...
##   ..$ VdnDe1_162036         : num [1:24] 1.49e-04 1.39e-04 8.97e-05 5.51e-05 5.29e-05 ...
##   ..$ PBdnRY1_2239          : num [1:24] 1.05e-04 7.82e-05 5.69e-05 1.17e-04 4.48e-05 ...
##   ..$ CLCdnPW22_15452       : num [1:24] 1.09e-04 1.01e-04 9.85e-05 6.06e-05 5.82e-05 ...
##   ..$ CLCdnPW50_2470        : num [1:24] NA NA NA NA NA NA NA NA NA NA ...
##   ..$ VdnPW3_37011          : num [1:24] 1.47e-04 7.78e-05 3.78e-05 1.16e-04 1.34e-04 ...
##   ..$ SdnPW1_46             : num [1:24] 6.09e-05 2.82e-05 2.74e-05 3.37e-05 1.30e-04 ...
##   ..$ CLCdnDe7_34340        : num [1:24] NA NA NA NA NA NA NA NA NA NA ...
##   ..$ CLCdnDe11_2286        : num [1:24] NA NA NA NA 0.000494 ...
##   ..$ VdnDe2_53978          : num [1:24] 0.000412 0.000536 0.000297 0.000548 0.000877 ...
##   ..$ CLCdnDe7_6144         : num [1:24] 0.000359 0.00025 0.000404 0.000199 0.000191 ...
##   ..$ CLCdnRY10_6741        : num [1:24] NA 1.34e-04 1.95e-04 NA 7.69e-05 ...
##   ..$ PBdnRY1_5389          : num [1:24] 0.00069 0.000495 0.000537 0.000626 0.001001 ...
##   ..$ VdnDe2_29290          : num [1:24] 0.000317 0.000134 0.000156 0.000192 0.000276 ...
##   ..$ VdnDe4_115100         : num [1:24] 0.000662 0.000673 0.000568 0.000628 0.001173 ...
##   ..$ Sotub11g025210.1.1    : num [1:24] 0.000272 0.000324 0.00021 0.000258 0.000289 ...
##   ..$ VdnDe1_39992          : num [1:24] 0.000529 0.000369 0.000417 0.00044 0.000282 ...
##   ..$ VdnPW4_15664          : num [1:24] NA NA NA NA NA NA NA NA 5.14e-05 6.16e-05 ...
##  $ metabolomics:'data.frame':    24 obs. of  22 variables:
##   ..$ Glukose : num [1:24] 1.41 0.91 1.16 1.56 2.07 2.56 3.18 5.01 1.01 6.43 ...
##   ..$ Fructose: num [1:24] 1.9 1.68 2.05 3.25 2.8 3.25 3.61 2.95 1.51 8.2 ...
##   ..$ Sucrose : num [1:24] 4.12 3.51 3.07 4.86 2.96 ...
##   ..$ Starch  : num [1:24] 5.47 4.65 5.44 6.2 2.03 ...
##   ..$ Asp     : num [1:24] 637 533 763 603 972 ...
##   ..$ Glu     : num [1:24] 3338 1824 3326 2785 2328 ...
##   ..$ Asn     : num [1:24] 273 249 496 404 1117 ...
##   ..$ Ser     : num [1:24] 520 303 364 439 289 ...
##   ..$ Gln     : num [1:24] 400 281 503 602 284 ...
##   ..$ Gly     : num [1:24] 81.8 59 69.8 81.9 89.1 65.4 94.9 97.6 45.7 71.3 ...
##   ..$ His     : num [1:24] 18.8 17.1 31.4 22.7 93.3 39.6 41.3 48.4 29.8 26.5 ...
##   ..$ Arg     : num [1:24] 31.9 22.3 56.2 46 124.1 ...
##   ..$ Thr     : num [1:24] 288 160 245 264 204 ...
##   ..$ Ala     : num [1:24] 1215 370 736 990 367 ...
##   ..$ Pro     : num [1:24] 104.9 52.8 83.5 122.3 120.4 ...
##   ..$ Tyr     : num [1:24] 27.5 24.1 37 32.9 114.6 ...
##   ..$ Val     : num [1:24] 64.9 35.6 73.9 73 86.3 ...
##   ..$ Met     : num [1:24] 7.4 4.5 6.2 7.9 7.6 6.6 6.7 6.2 3.5 2.6 ...
##   ..$ Ile     : num [1:24] 82.9 49.6 93 80.1 190.8 ...
##   ..$ Lys     : num [1:24] 31.1 23 44.5 28.1 39.9 27.4 40.8 47 35.9 25.6 ...
##   ..$ Leu     : num [1:24] 26.8 13.7 34.3 26.4 34.3 22.4 23.3 35.5 43.4 32.2 ...
##   ..$ Phe     : num [1:24] 123.6 76.3 122.9 121.7 169.9 ...
##  $ hormonomics :'data.frame':    24 obs. of  12 variables:
##   ..$ IAA       : num [1:24] 49.8 35.7 37.9 51.4 86.6 ...
##   ..$ oxIAA     : num [1:24] 51.8 52.8 53 68.3 49.2 ...
##   ..$ IAA.Asp   : num [1:24] 2.04 1.55 2.53 2.02 1.75 ...
##   ..$ ABA       : num [1:24] 47.4 51.3 42.8 48.6 43 ...
##   ..$ PA        : num [1:24] 84 124.6 99.6 114.1 135 ...
##   ..$ DPA       : num [1:24] 94.8 108 88.1 122 265.4 ...
##   ..$ SA        : num [1:24] 524 591 448 948 1113 ...
##   ..$ JA        : num [1:24] 7.07 5.58 4.21 7.67 9.58 ...
##   ..$ JA.Ile    : num [1:24] 0.663 0.325 0.551 0.34 0.845 ...
##   ..$ X9.10.dhJA: num [1:24] 5.68 8.21 8.92 5.7 7.81 ...
##   ..$ X12.OH.JA : num [1:24] 70.3 53.2 42.3 87.2 185.4 ...
##   ..$ cisOPDA   : num [1:24] 518 440 468 369 1382 ...
##  $ qPCR        :'data.frame':    24 obs. of  14 variables:
##   ..$ RbohA  : num [1:24] 0.687 0.517 0.877 0.686 0.32 ...
##   ..$ SnRK2  : num [1:24] 1.811 1.578 2.114 1.757 0.773 ...
##   ..$ ACO2   : num [1:24] 0.1086 0.1314 0.1291 0.1591 0.0584 ...
##   ..$ HSP70  : num [1:24] 0.798 0.774 1.127 0.856 0.672 ...
##   ..$ PR1b   : num [1:24] 0.2356 0.1818 0.1353 0.1406 0.0485 ...
##   ..$ RD29B  : num [1:24] 0.017 0.017 0.017 0.022 0.113 ...
##   ..$ X13.LOX: num [1:24] 0.791 0.612 0.771 0.743 0.428 ...
##   ..$ P5CS   : num [1:24] 1.959 1.586 2.249 1.535 0.995 ...
##   ..$ ERF1   : num [1:24] 0.689 0.503 0.563 0.793 0.517 ...
##   ..$ CAT1   : num [1:24] 0.949 0.814 0.699 0.632 0.222 ...
##   ..$ CO     : num [1:24] 2.189 1.879 2.548 1.791 0.911 ...
##   ..$ SWEET  : num [1:24] 1.301 1.073 1.391 1.068 0.449 ...
##   ..$ SP6A   : num [1:24] 0.1224 0.1224 0.1224 0.1519 0.0659 ...
##   ..$ M0ZJG3 : num [1:24] 0.9759 0.7177 0.0848 1.0922 1.0931 ...
##  $ phenomics   :'data.frame':    24 obs. of  8 variables:
##   ..$ FW_SH_Total      : num [1:24] 52.2 91.7 60.1 93.3 92.4 ...
##   ..$ DW_SH_Total      : num [1:24] 2.78 4.7 3.21 4.68 5.78 ...
##   ..$ TOP.AREA         : num [1:24] 96175 95186 96925 103507 119552 ...
##   ..$ COMPACTNESS      : num [1:24] 0.742 0.728 0.724 0.721 0.781 ...
##   ..$ Fv.Fm            : num [1:24] 0.719 0.717 0.717 0.722 0.711 ...
##   ..$ qL               : num [1:24] 0.798 0.789 0.768 0.796 0.776 ...
##   ..$ deltaT           : num [1:24] 0.7062 0.0231 0.7603 0 0.7342 ...
##   ..$ WATER.CONSUMPTION: int [1:24] 288 114 119 130 285 211 199 182 205 210 ...

names(X)

## [1] "status"       "proteomics"   "metabolomics" "hormonomics" 
## [5] "qPCR"         "phenomics"

Check if sample names in all datasets match.

OK <- TRUE
for (i in 2:length(X)) {
    print(ok <- all(names(X[[1]]) == rownames(X[[i]])))
    OK <- OK & ok
}

## [1] TRUE
## [1] TRUE
## [1] TRUE
## [1] TRUE
## [1] TRUE

Sample names in datasets match.

2.5.1 Check missing values

Pairs of Variables with many missing values sometimes produce missing correlations. Most of the times they will be ignored in the analysis but might produce problems at some point.

This checks if any such pairs exist:

for (i in 2:length(X)) {
    catln(names(X[i]))
    bad_columns <- find_NA_pairs(X[[i]])
    if (!is.null(bad_columns)) {
        cat("   Columns with mutual NA values:\n")
        cat(paste("     ", bad_columns, collapse = ",\n"), sep = "\n")
    }
}

## proteomics 
##    Columns with mutual NA values:
##       VdnPW1_81435,
##       CLCdnDe13_5738,
##       CLCdnDe10_24187,
##       CLCdnDe7_34340,
##       CLCdnDe2_96512,
##       PBdnRY1_11085,
##       CLCdnPW50_2470,
##       CLCdnDe11_2286,
##       VdnPW4_15664
## metabolomics 
## hormonomics 
## qPCR 
## phenomics

Listed vector of names (if not NULL) can be used to delete the columns with lots of missing values.

Put data into safe object DATA.

DATA <- X

We will also need the names of treatment groups.

groups <- unique(pdata$Treatment)
groups

## [1] "H"  "HD"

3 Multiomics data integration with DIABLO

CCDATA <- DATA
names(CCDATA)

## [1] "status"       "proteomics"   "metabolomics" "hormonomics" 
## [5] "qPCR"         "phenomics"

out <- ""
out <- paste(out, knit_child("035-DIABLO-2.Rmd", quiet = TRUE))

## Warning in eng_r(options): Failed to tidy R code in chunk 'plotArrow'. Reason:
## Error in base::parse(text = code, keep.source = FALSE) : 
##   <text>:3:1: unexpected symbol
## 2: title = paste ( groups , collapse = ", " )
## 3: invisible
##    ^

cat(out)

Child: 035-DIABLO-2.Rmd ## DIABLO proteomics, metabolomics, hormonomics, qPCR, phenomics ## DIABLO

DIABLO from mixOmics enables integration of more than two datasets.

3.1 Organize data

Thre datasets are organized as a list of matrices with same samples as rows and variables in columns.

data <- CCDATA[-1]
str(data)

## List of 5
##  $ proteomics  :'data.frame':    24 obs. of  36 variables:
##   ..$ VdnDe2_86784          : num [1:24] 0.000281 0.000522 0.00038 0.000467 0.000374 ...
##   ..$ PBdnRY1_7025          : num [1:24] 0.000478 0.000443 0.000603 0.000529 0.000407 ...
##   ..$ VdnPW5_1541           : num [1:24] 0.00131 0.00218 0.00212 0.00152 0.00104 ...
##   ..$ CLCdnPW10_11437       : num [1:24] 0.000612 0.000455 0.000441 0.000543 0.000782 ...
##   ..$ CLCdnDe2_96512        : num [1:24] NA NA NA NA NA ...
##   ..$ VdnPW1_81435          : num [1:24] 2.87e-05 2.66e-05 2.58e-05 3.18e-05 NA NA NA NA NA NA ...
##   ..$ CLCdnDe6_4953         : num [1:24] 0.000569 0.000447 0.000395 0.000485 0.000512 ...
##   ..$ CLCdnDe13_5738        : num [1:24] NA NA NA NA 0.000149 ...
##   ..$ CLCdnPW12_1472        : num [1:24] 0.000542 0.000503 0.000466 0.000545 0.000483 ...
##   ..$ CLCdnPW49_8560        : num [1:24] 3.97e-05 3.69e-05 NA NA 8.45e-05 ...
##   ..$ VdnPW4_18801          : num [1:24] 0.000276 0.000192 0.000218 0.000268 0.000404 ...
##   ..$ SdnRY1_10712          : num [1:24] 0.00037 0.000386 0.000375 0.000359 0.000295 ...
##   ..$ CLCdnDe2_311          : num [1:24] 0.000262 0.000122 0.000157 0.000194 0.000186 ...
##   ..$ CLCdnDe13_63179       : num [1:24] 9.19e-05 8.53e-05 4.14e-05 2.55e-04 1.96e-04 ...
##   ..$ CLCdnDe10_24187       : num [1:24] NA NA NA NA 0.00013 ...
##   ..$ TDe_comp110944_c2_seq4: num [1:24] 5.15e-05 4.78e-05 NA 5.71e-05 1.64e-04 ...
##   ..$ PBdnRY1_11085         : num [1:24] NA NA NA NA NA NA NA NA 6.20e-05 7.43e-05 ...
##   ..$ CLCdnPW48_15164       : num [1:24] 0.000201 0.00028 0.000333 0.000186 0.000357 ...
##   ..$ CLCdnPW31_956         : num [1:24] 0.000136 0.00019 0.000184 0.000251 0.000217 ...
##   ..$ VdnDe1_162036         : num [1:24] 1.49e-04 1.39e-04 8.97e-05 5.51e-05 5.29e-05 ...
##   ..$ PBdnRY1_2239          : num [1:24] 1.05e-04 7.82e-05 5.69e-05 1.17e-04 4.48e-05 ...
##   ..$ CLCdnPW22_15452       : num [1:24] 1.09e-04 1.01e-04 9.85e-05 6.06e-05 5.82e-05 ...
##   ..$ CLCdnPW50_2470        : num [1:24] NA NA NA NA NA NA NA NA NA NA ...
##   ..$ VdnPW3_37011          : num [1:24] 1.47e-04 7.78e-05 3.78e-05 1.16e-04 1.34e-04 ...
##   ..$ SdnPW1_46             : num [1:24] 6.09e-05 2.82e-05 2.74e-05 3.37e-05 1.30e-04 ...
##   ..$ CLCdnDe7_34340        : num [1:24] NA NA NA NA NA NA NA NA NA NA ...
##   ..$ CLCdnDe11_2286        : num [1:24] NA NA NA NA 0.000494 ...
##   ..$ VdnDe2_53978          : num [1:24] 0.000412 0.000536 0.000297 0.000548 0.000877 ...
##   ..$ CLCdnDe7_6144         : num [1:24] 0.000359 0.00025 0.000404 0.000199 0.000191 ...
##   ..$ CLCdnRY10_6741        : num [1:24] NA 1.34e-04 1.95e-04 NA 7.69e-05 ...
##   ..$ PBdnRY1_5389          : num [1:24] 0.00069 0.000495 0.000537 0.000626 0.001001 ...
##   ..$ VdnDe2_29290          : num [1:24] 0.000317 0.000134 0.000156 0.000192 0.000276 ...
##   ..$ VdnDe4_115100         : num [1:24] 0.000662 0.000673 0.000568 0.000628 0.001173 ...
##   ..$ Sotub11g025210.1.1    : num [1:24] 0.000272 0.000324 0.00021 0.000258 0.000289 ...
##   ..$ VdnDe1_39992          : num [1:24] 0.000529 0.000369 0.000417 0.00044 0.000282 ...
##   ..$ VdnPW4_15664          : num [1:24] NA NA NA NA NA NA NA NA 5.14e-05 6.16e-05 ...
##  $ metabolomics:'data.frame':    24 obs. of  22 variables:
##   ..$ Glukose : num [1:24] 1.41 0.91 1.16 1.56 2.07 2.56 3.18 5.01 1.01 6.43 ...
##   ..$ Fructose: num [1:24] 1.9 1.68 2.05 3.25 2.8 3.25 3.61 2.95 1.51 8.2 ...
##   ..$ Sucrose : num [1:24] 4.12 3.51 3.07 4.86 2.96 ...
##   ..$ Starch  : num [1:24] 5.47 4.65 5.44 6.2 2.03 ...
##   ..$ Asp     : num [1:24] 637 533 763 603 972 ...
##   ..$ Glu     : num [1:24] 3338 1824 3326 2785 2328 ...
##   ..$ Asn     : num [1:24] 273 249 496 404 1117 ...
##   ..$ Ser     : num [1:24] 520 303 364 439 289 ...
##   ..$ Gln     : num [1:24] 400 281 503 602 284 ...
##   ..$ Gly     : num [1:24] 81.8 59 69.8 81.9 89.1 65.4 94.9 97.6 45.7 71.3 ...
##   ..$ His     : num [1:24] 18.8 17.1 31.4 22.7 93.3 39.6 41.3 48.4 29.8 26.5 ...
##   ..$ Arg     : num [1:24] 31.9 22.3 56.2 46 124.1 ...
##   ..$ Thr     : num [1:24] 288 160 245 264 204 ...
##   ..$ Ala     : num [1:24] 1215 370 736 990 367 ...
##   ..$ Pro     : num [1:24] 104.9 52.8 83.5 122.3 120.4 ...
##   ..$ Tyr     : num [1:24] 27.5 24.1 37 32.9 114.6 ...
##   ..$ Val     : num [1:24] 64.9 35.6 73.9 73 86.3 ...
##   ..$ Met     : num [1:24] 7.4 4.5 6.2 7.9 7.6 6.6 6.7 6.2 3.5 2.6 ...
##   ..$ Ile     : num [1:24] 82.9 49.6 93 80.1 190.8 ...
##   ..$ Lys     : num [1:24] 31.1 23 44.5 28.1 39.9 27.4 40.8 47 35.9 25.6 ...
##   ..$ Leu     : num [1:24] 26.8 13.7 34.3 26.4 34.3 22.4 23.3 35.5 43.4 32.2 ...
##   ..$ Phe     : num [1:24] 123.6 76.3 122.9 121.7 169.9 ...
##  $ hormonomics :'data.frame':    24 obs. of  12 variables:
##   ..$ IAA       : num [1:24] 49.8 35.7 37.9 51.4 86.6 ...
##   ..$ oxIAA     : num [1:24] 51.8 52.8 53 68.3 49.2 ...
##   ..$ IAA.Asp   : num [1:24] 2.04 1.55 2.53 2.02 1.75 ...
##   ..$ ABA       : num [1:24] 47.4 51.3 42.8 48.6 43 ...
##   ..$ PA        : num [1:24] 84 124.6 99.6 114.1 135 ...
##   ..$ DPA       : num [1:24] 94.8 108 88.1 122 265.4 ...
##   ..$ SA        : num [1:24] 524 591 448 948 1113 ...
##   ..$ JA        : num [1:24] 7.07 5.58 4.21 7.67 9.58 ...
##   ..$ JA.Ile    : num [1:24] 0.663 0.325 0.551 0.34 0.845 ...
##   ..$ X9.10.dhJA: num [1:24] 5.68 8.21 8.92 5.7 7.81 ...
##   ..$ X12.OH.JA : num [1:24] 70.3 53.2 42.3 87.2 185.4 ...
##   ..$ cisOPDA   : num [1:24] 518 440 468 369 1382 ...
##  $ qPCR        :'data.frame':    24 obs. of  14 variables:
##   ..$ RbohA  : num [1:24] 0.687 0.517 0.877 0.686 0.32 ...
##   ..$ SnRK2  : num [1:24] 1.811 1.578 2.114 1.757 0.773 ...
##   ..$ ACO2   : num [1:24] 0.1086 0.1314 0.1291 0.1591 0.0584 ...
##   ..$ HSP70  : num [1:24] 0.798 0.774 1.127 0.856 0.672 ...
##   ..$ PR1b   : num [1:24] 0.2356 0.1818 0.1353 0.1406 0.0485 ...
##   ..$ RD29B  : num [1:24] 0.017 0.017 0.017 0.022 0.113 ...
##   ..$ X13.LOX: num [1:24] 0.791 0.612 0.771 0.743 0.428 ...
##   ..$ P5CS   : num [1:24] 1.959 1.586 2.249 1.535 0.995 ...
##   ..$ ERF1   : num [1:24] 0.689 0.503 0.563 0.793 0.517 ...
##   ..$ CAT1   : num [1:24] 0.949 0.814 0.699 0.632 0.222 ...
##   ..$ CO     : num [1:24] 2.189 1.879 2.548 1.791 0.911 ...
##   ..$ SWEET  : num [1:24] 1.301 1.073 1.391 1.068 0.449 ...
##   ..$ SP6A   : num [1:24] 0.1224 0.1224 0.1224 0.1519 0.0659 ...
##   ..$ M0ZJG3 : num [1:24] 0.9759 0.7177 0.0848 1.0922 1.0931 ...
##  $ phenomics   :'data.frame':    24 obs. of  8 variables:
##   ..$ FW_SH_Total      : num [1:24] 52.2 91.7 60.1 93.3 92.4 ...
##   ..$ DW_SH_Total      : num [1:24] 2.78 4.7 3.21 4.68 5.78 ...
##   ..$ TOP.AREA         : num [1:24] 96175 95186 96925 103507 119552 ...
##   ..$ COMPACTNESS      : num [1:24] 0.742 0.728 0.724 0.721 0.781 ...
##   ..$ Fv.Fm            : num [1:24] 0.719 0.717 0.717 0.722 0.711 ...
##   ..$ qL               : num [1:24] 0.798 0.789 0.768 0.796 0.776 ...
##   ..$ deltaT           : num [1:24] 0.7062 0.0231 0.7603 0 0.7342 ...
##   ..$ WATER.CONSUMPTION: int [1:24] 288 114 119 130 285 211 199 182 205 210 ...

length(data)

## [1] 5

In addition, outcome, phenotypic state or in our case treatment can also be determined. We have combination of two treatments and four time points.

state <- factor(CCDATA[[1]])
table(state)

## state
##  H01  H07  H08  H14 HD08 HD14 
##    4    4    4    4    4    4

str(state)

##  Factor w/ 6 levels "H01","H07","H08",..: 1 1 1 1 2 2 2 2 3 3 ...
##  - attr(*, "names")= chr [1:24] "H_S1_1" "H_S1_2" "H_S1_3" "H_S1_4" ...

3.2 Initial analysis

Note: Additional insights can be found at http://mixomics.org/mixdiablo/diablo-tcga-case-study/.

3.2.1 Pairwise PLS comparisons

list.keepX = c(25, 25)  # select arbitrary values of features to keep
list.keepY = c(25, 25)

par(mfrow = c(2, 2))
pairs <- combn(1:length(names(data)), 2)
nms <- names(data)
outn <- ""
j <- 1
ncomp <- length(data)

cols <- c("orange1", "brown1", "lightgreen", "lightblue", "pink")[1:length(data)]
pchs <- c(16, 17, 15, 18, 20)[1:length(data)]
j <- 1
cutoff <- 0.5
for (j in 1:ncol(pairs)) {
    pair <- pairs[, j]
    pair
    X <- CCDATA[[pair[1] + 1]]
    Y <- CCDATA[[pair[2] + 1]]
    list.keepX <- rep(min(ncol(X), 25), ncomp)
    list.keepY <- rep(min(ncol(Y), 25), ncomp)
    x <- spls(X, Y, ncomp = ncomp, keepX = list.keepX, keepY = list.keepY)
    assign(paste0("spls", j), x)
    cat("\n", paste(nms[pair]), "\n")
    cat("Results in:", paste0("spls", j), "\n")
    cat("Correlation between pls variates:\n")
    print(round(cor(x$variates$X, x$variates$Y), 5))
    # 
    plotVar(x, cutoff = cutoff, title = paste(nms[pair], collapse = ", "), 
        legend = c(nms[pair][1], nms[pair][2]), var.names = FALSE, style = "graphics", 
        pch = pchs[pair], cex = c(2, 2), col = cols[pair])
}

## 
##  proteomics metabolomics 
## Results in: spls1 
## Correlation between pls variates:
##          comp1    comp2    comp3    comp4    comp5
## comp1  0.83322 -0.00024  0.00000 -0.00004 -0.00212
## comp2 -0.26131  0.79508  0.00000  0.00003  0.00719
## comp3 -0.00278 -0.09338  0.83087  0.00011 -0.00179
## comp4 -0.11951  0.00605 -0.09477  0.74209  0.00026
## comp5  0.01624 -0.03563 -0.20345  0.26945  0.64341

## 
##  proteomics hormonomics 
## Results in: spls2 
## Correlation between pls variates:
##          comp1    comp2   comp3    comp4    comp5
## comp1  0.89069 -0.00009 0.00118  0.00236  0.00827
## comp2 -0.12695  0.84463 0.00003 -0.00051 -0.00112
## comp3  0.11435 -0.09072 0.82335 -0.00007  0.00070
## comp4 -0.18493 -0.12829 0.09514  0.89464  0.00005
## comp5  0.13655 -0.25183 0.17718  0.06506  0.84667

## 
##  proteomics qPCR 
## Results in: spls3 
## Correlation between pls variates:
##          comp1    comp2    comp3    comp4    comp5
## comp1  0.77134  0.00277 -0.02179  0.01149  0.01837
## comp2  0.44964  0.72964 -0.00125  0.00132 -0.00204
## comp3 -0.22312 -0.38778  0.67382 -0.00004 -0.00125
## comp4  0.06602  0.02664  0.19081  0.81444 -0.00051
## comp5  0.20256  0.28200 -0.36973  0.21080  0.62801

## 
##  proteomics phenomics 
## Results in: spls4 
## Correlation between pls variates:
##          comp1    comp2    comp3    comp4    comp5
## comp1  0.84289 -0.00006  0.00220 -0.00341  0.02694
## comp2 -0.18976  0.88902  0.00010  0.00007 -0.00637
## comp3  0.20550  0.06472  0.83362  0.00002  0.00973
## comp4  0.35242 -0.13487 -0.00499  0.75882 -0.00157
## comp5 -0.18724  0.15666  0.07214 -0.37875  0.60641

Circle Correlation Plots for pairwise PLS models on ADAPT data. At most top 25 features for each dimension with correlation above 0.5, are displayed.

## 
##  metabolomics hormonomics 
## Results in: spls5 
## Correlation between pls variates:
##          comp1    comp2    comp3    comp4   comp5
## comp1  0.92070  0.00000  0.00000  0.00000 0.00000
## comp2 -0.22655  0.80614  0.00000  0.00000 0.00000
## comp3  0.03801 -0.09867  0.80378  0.00000 0.00000
## comp4  0.01947  0.06284  0.00154  0.75077 0.00000
## comp5  0.06147  0.15645 -0.29720 -0.20970 0.73359

## 
##  metabolomics qPCR 
## Results in: spls6 
## Correlation between pls variates:
##          comp1    comp2    comp3    comp4  comp5
## comp1  0.56114  0.00000  0.00000  0.00000 0.0000
## comp2  0.59449  0.71184  0.00000  0.00000 0.0000
## comp3  0.31180  0.35694  0.54464  0.00000 0.0000
## comp4  0.26483  0.31609  0.44040  0.52048 0.0000
## comp5 -0.13188 -0.17760 -0.24392 -0.28950 0.3771

## 
##  metabolomics phenomics 
## Results in: spls7 
## Correlation between pls variates:
##          comp1    comp2    comp3    comp4   comp5
## comp1  0.85036  0.00000  0.00000  0.00000 0.00000
## comp2  0.39864  0.79420  0.00000  0.00000 0.00000
## comp3 -0.08328 -0.34091  0.71280  0.00000 0.00000
## comp4  0.07018  0.03387 -0.06439  0.76130 0.00000
## comp5  0.14899  0.01122 -0.03306 -0.20873 0.83194

## 
##  hormonomics qPCR 
## Results in: spls8 
## Correlation between pls variates:
##          comp1    comp2    comp3   comp4   comp5
## comp1  0.67791  0.00000  0.00000 0.00000 0.00000
## comp2 -0.34141  0.48630  0.00000 0.00000 0.00000
## comp3 -0.20327  0.29962  0.41743 0.00000 0.00000
## comp4  0.13282 -0.09944 -0.09450 0.43011 0.00000
## comp5  0.30685 -0.39619 -0.42308 0.46567 0.50213

Circle Correlation Plots for pairwise PLS models on ADAPT data. At most top 25 features for each dimension with correlation above 0.5, are displayed.

## 
##  hormonomics phenomics 
## Results in: spls9 
## Correlation between pls variates:
##          comp1    comp2    comp3   comp4   comp5
## comp1  0.86662  0.00000  0.00000 0.00000 0.00000
## comp2 -0.16299  0.77552  0.00000 0.00000 0.00000
## comp3 -0.03923  0.05898  0.65368 0.00000 0.00000
## comp4  0.01873 -0.22032 -0.37616 0.67806 0.00000
## comp5  0.18418 -0.13481 -0.07150 0.12244 0.44097

## 
##  qPCR phenomics 
## Results in: spls10 
## Correlation between pls variates:
##          comp1    comp2    comp3   comp4   comp5
## comp1  0.69225  0.00000  0.00000 0.00000 0.00000
## comp2  0.36644  0.53084  0.00000 0.00000 0.00000
## comp3 -0.15153 -0.15885  0.51075 0.00000 0.00000
## comp4  0.18052  0.28591 -0.13513 0.35046 0.00000
## comp5  0.05498  0.09618 -0.14546 0.15602 0.38489

Circle Correlation Plots for pairwise PLS models on ADAPT data. At most top 25 features for each dimension with correlation above 0.5, are displayed.

3.2.2 Initial DIABLO model

Following the suggestion in the source, we will use design matrices with small values. This is supposed to keep low classification error rate.

entry <- .entry
design = matrix(entry, ncol = length(data), nrow = length(data), dimnames = list(names(data), 
    names(data)))
diag(design) = 0  # set diagonal to 0s

design

##              proteomics metabolomics hormonomics qPCR phenomics
## proteomics          0.0          0.5         0.5  0.5       0.5
## metabolomics        0.5          0.0         0.5  0.5       0.5
## hormonomics         0.5          0.5         0.0  0.5       0.5
## qPCR                0.5          0.5         0.5  0.0       0.5
## phenomics           0.5          0.5         0.5  0.5       0.0

With a design in place, the initial DIABLO model can be generated. An arbitrarily high number of components (ncomp = 5) will be used.

# form basic DIABLO model
Y <- state
basic.diablo.model = block.splsda(X = data, Y = Y, ncomp = 5, design = design)

## Design matrix has changed to include Y; each block will be
##             linked to Y.

######################################## test eval
knitr::opts_chunk$set(eval = !FALSE)

3.2.3 Tuning the number of components

Details of tuning process can be found in the http://mixomics.org/mixdiablo/diablo-tcga-case-study/.

The process can be computer time consuming and was performed separately.

From previous tuning sessions one can conclude, that the classification rate stays roughly unchanged after two to four components, so we will set the number of components to the number of data sets:

ncomp <- length(data)
ncomp

## [1] 5

3.2.4 Tuning the number of features

We choose the optimal number of variables to select in each data set using the tune.block.splsda() function, for a grid of keepX values for each type of omics. Note that the function has been set to favour a relatively small signature while allowing us to obtain a sufficient number of variables for downstream validation and/or interpretation. See ?tune.block.splsda.

Previous analyses suggest the following list:

$metabolomics
[1] 10 10 5
$hormonomics
[1] 5 5 10
$qPCR
[1] 10 10 5

Number of kept variates is limited with the minimal number of variates in the datasets. We have decided to keep at most 10 variates for each component.

min_variates <- min(sapply(data, ncol), 10)
if (FALSE) keepX <- list(metabolomics = rep(10, ncomp), hormonomics = rep(10, 
    ncomp), qPCR = rep(10, ncomp))
list.keepX = list()
for (i in 1:length(data)) list.keepX[[i]] <- rep(min_variates, ncomp)
names(list.keepX) <- names(data)
list.keepX

## $proteomics
## [1] 8 8 8 8 8
## 
## $metabolomics
## [1] 8 8 8 8 8
## 
## $hormonomics
## [1] 8 8 8 8 8
## 
## $qPCR
## [1] 8 8 8 8 8
## 
## $phenomics
## [1] 8 8 8 8 8

3.3 Final model

The final DIABLO model is run as:

# set the optimised DIABLO model
final.diablo.model = block.splsda(X = data, Y = as.factor(state), ncomp = ncomp, 
    keepX = list.keepX, design = design)

## Design matrix has changed to include Y; each block will be
##             linked to Y.

The selected variables can be extracted with the function selectVar(), for example in each block, as seen below. Note that the stability of selected variables can be extracted from the output of the perf() function.

# the features selected from components
for (comp in 1:ncomp) {
    cat("\nComponent ", comp, ":\n")
    for (i in 1:length(data)) {
        cat(names(data)[i], "\n")
        print(selectVar(final.diablo.model, comp = comp)[[i]]$name)
    }
}

## 
## Component  1 :
## proteomics 
## [1] "PBdnRY1_2239"    "PBdnRY1_5389"    "PBdnRY1_7025"   
## [4] "CLCdnDe6_4953"   "CLCdnPW49_8560"  "VdnPW3_37011"   
## [7] "SdnPW1_46"       "CLCdnPW10_11437"
## metabolomics 
## [1] "Glukose"  "Fructose" "Ile"      "Starch"   "Val"      "Gly"     
## [7] "His"      "Tyr"     
## hormonomics 
## [1] "SA"         "ABA"        "X12.OH.JA"  "DPA"        "JA.Ile"    
## [6] "PA"         "cisOPDA"    "X9.10.dhJA"
## qPCR 
## [1] "RD29B"   "X13.LOX" "SP6A"    "M0ZJG3"  "RbohA"   "HSP70"  
## [7] "CAT1"    "PR1b"   
## phenomics 
## [1] "Fv.Fm"             "DW_SH_Total"       "deltaT"           
## [4] "qL"                "FW_SH_Total"       "TOP.AREA"         
## [7] "COMPACTNESS"       "WATER.CONSUMPTION"
## 
## Component  2 :
## proteomics 
## [1] "VdnDe2_53978"       "CLCdnPW48_15164"    "CLCdnPW12_1472"    
## [4] "Sotub11g025210.1.1" "CLCdnRY10_6741"     "CLCdnPW31_956"     
## [7] "VdnPW4_15664"       "CLCdnPW22_15452"   
## metabolomics 
## [1] "Ala" "Gln" "Met" "Asn" "Thr" "Leu" "Gly" "Arg"
## hormonomics 
## [1] "oxIAA"      "JA"         "X9.10.dhJA" "IAA.Asp"    "ABA"       
## [6] "IAA"        "X12.OH.JA"  "DPA"       
## qPCR 
## [1] "ACO2"  "SnRK2" "ERF1"  "PR1b"  "RbohA" "HSP70" "CAT1"  "P5CS" 
## phenomics 
## [1] "TOP.AREA"          "WATER.CONSUMPTION" "COMPACTNESS"      
## [4] "Fv.Fm"             "FW_SH_Total"       "deltaT"           
## [7] "DW_SH_Total"       "qL"               
## 
## Component  3 :
## proteomics 
## [1] "CLCdnDe7_6144"   "SdnRY1_10712"    "VdnDe1_39992"   
## [4] "SdnPW1_46"       "CLCdnPW49_8560"  "CLCdnPW31_956"  
## [7] "VdnPW4_18801"    "CLCdnDe13_63179"
## metabolomics 
## [1] "Asp"     "Leu"     "Lys"     "Ser"     "Pro"     "His"    
## [7] "Sucrose" "Arg"    
## hormonomics 
## [1] "X12.OH.JA"  "X9.10.dhJA" "cisOPDA"    "SA"         "PA"        
## [6] "JA"         "JA.Ile"     "DPA"       
## qPCR 
## [1] "RbohA"   "PR1b"    "RD29B"   "ERF1"    "P5CS"    "SP6A"   
## [7] "ACO2"    "X13.LOX"
## phenomics 
## [1] "COMPACTNESS"       "qL"                "deltaT"           
## [4] "WATER.CONSUMPTION" "DW_SH_Total"       "FW_SH_Total"      
## [7] "Fv.Fm"             "TOP.AREA"         
## 
## Component  4 :
## proteomics 
## [1] "CLCdnPW10_11437" "VdnDe4_115100"   "VdnPW4_18801"   
## [4] "VdnPW5_1541"     "CLCdnPW22_15452" "VdnDe2_86784"   
## [7] "CLCdnDe10_24187" "CLCdnPW12_1472" 
## metabolomics 
## [1] "Arg"     "Phe"     "His"     "Asn"     "Asp"     "Met"    
## [7] "Sucrose" "Tyr"    
## hormonomics 
## [1] "IAA"        "IAA.Asp"    "JA.Ile"     "PA"         "X9.10.dhJA"
## [6] "JA"         "X12.OH.JA"  "ABA"       
## qPCR 
## [1] "P5CS"    "ERF1"    "M0ZJG3"  "PR1b"    "SP6A"    "X13.LOX"
## [7] "RD29B"   "RbohA"  
## phenomics 
## [1] "WATER.CONSUMPTION" "Fv.Fm"             "FW_SH_Total"      
## [4] "DW_SH_Total"       "qL"                "deltaT"           
## [7] "COMPACTNESS"       "TOP.AREA"         
## 
## Component  5 :
## proteomics 
## [1] "VdnDe2_86784"           "TDe_comp110944_c2_seq4"
## [3] "VdnPW5_1541"            "Sotub11g025210.1.1"    
## [5] "CLCdnDe2_311"           "PBdnRY1_2239"          
## [7] "VdnDe1_162036"          "CLCdnPW22_15452"       
## metabolomics 
## [1] "Glu" "Asn" "Leu" "Ser" "His" "Pro" "Tyr" "Thr"
## hormonomics 
## [1] "PA"        "JA.Ile"    "cisOPDA"   "IAA"       "X12.OH.JA"
## [6] "SA"        "JA"        "ABA"      
## qPCR 
## [1] "PR1b"    "SP6A"    "X13.LOX" "RD29B"   "CAT1"    "ERF1"   
## [7] "SWEET"   "M0ZJG3" 
## phenomics 
## [1] "qL"                "DW_SH_Total"       "COMPACTNESS"      
## [4] "WATER.CONSUMPTION" "deltaT"            "FW_SH_Total"      
## [7] "TOP.AREA"          "Fv.Fm"

3.4 Sample plots

plotDIABLO() is a diagnostic plot to check whether the correlation between components from each data set has been maximised as specified in the design matrix. We specify which dimension to be assessed with the ncomp argument.

for (comp in 1:ncomp) {
    plotDiablo(final.diablo.model, ncomp = comp)
    title(paste("Component", comp), adj = 0.1, line = -1, outer = TRUE)
}

The sample plot with the plotIndiv() function projects each sample into the space spanned by the components of each block. Clustering of the samples can be assessed with this plot.

plind <- plotIndiv(final.diablo.model, ind.names = FALSE, legend = TRUE, 
    title = "DIABLO Sample Plots", guide = "none", ellipse = TRUE)

## Warning: It is deprecated to specify `guide = FALSE` to remove a
## guide. Please use `guide = "none"` instead.

In the arrow plot below, the start of the arrow indicates the centroid between all data sets for a given sample and the tips of the arrows indicate the location of that sample in each block. Such graphics highlight the agreement between all data sets at the sample level. While somewhat difficult to interpret, even qualitatively, this arrow plot shows proximities of C01 and H01 (both day 1), C07 and C08, and H07 and H08 ( both a day apart). While C samples are in forth quadrant ( D1 < 0, D2 > 0), H samples have ( D1 < 0, D2 < 0) except H14 that is separated on the positive part of Dimension 1.

plotArrow(final.diablo.model, ind.names = FALSE, legend = TRUE,
          title = paste(groups,collapse=", ")

          )

3.5 Variable plots

Several graphical outputs are available to visualise and mine the associations between the selected variables.

The best starting point to evaluate the correlation structure between variables is with the correlation circle plot. A majority of the qPCR variables are positively correlated only with the first component. The hormonomics and metabolomics variables seem to separate along first two dimensions. These first two components correlate highly with the selected variables from the all three dataset. From this, the correlation of each selected feature from all three datasets can be evaluated based on their proximity.

# if(length(data)==3) pick <- 1:3 else pick <- c(4,1:3) cols <-
# c('orange1', 'brown1', 'lightgreen','lightblue')[pick] pchs <- c(16,
# 17, 15, 18)[pick]
cols <- c("orange1", "brown1", "lightgreen", "lightblue", "pink")[1:length(data)]
pchs <- c(16, 17, 15, 18, 20)[1:length(data)]
plotVar(final.diablo.model, var.names = FALSE, style = "graphics", legend = TRUE, 
    pch = pchs, cex = rep(2, length(data)), col = cols)

The circos plot is exclusive to integrative frameworks and represents the correlations between variables of different types, represented on the side quadrants. It seems that the qPCR variables are almost entirely negatively correlated with the other two dataframes. Just few from metabolomics and hormonomics are positively correlated. Note that these correlations are above a value of 0.7 (cutoff = 0.7). All interpretations made above are only relevant for features with very strong correlations.

Plot variables

# plotVar(res, cutoff=0.5, legend = TRUE, overlap=!FALSE,
# style='graphics') plotVar(res, cutoff=0.5, legend = TRUE,
# overlap=FALSE, style='graphics')
plotVar(final.diablo.model, cutoff = 0.5, legend = FALSE, comp = c(1, 2), 
    overlap = FALSE, col = cols, cex = rep(4, length(data)))

plotVar(final.diablo.model, cutoff = 0.5, legend = FALSE, comp = c(2, 3), 
    overlap = FALSE, col = cols, cex = rep(4, length(data)))

circosPlot(final.diablo.model, cutoff = 0.7, line = TRUE, color.blocks = cols, 
    color.cor = c(3, 2), size.labels = 1, xpd = TRUE)

3.6 Relevance networks

Another visualisation of the correlations between the different types of variables is the relevance network, which is also built on the similarity matrix (as is the circos plot). Colour represent variable type.

blocks <- combn(length(data), 2)
j <- 1
cutoff <- 0.7
out35a <- ""
for (j in 1:ncol(blocks)) {
    out35a <- paste(out35a, knit_child("035a-DIABLO-network.Rmd", quiet = !TRUE))
    
    if (interactive()) 
        readline()
}
cat(out35a)

{r echo=FALSE}gsub("\\.","_",as.character(cutoff)############################################### ## ## ## (c) Andrej Blejec (andrej.blejec@nib.si) ## ## ## ###############################################

3.6.1 proteomics and metabolomics

nfn <- paste0("../reports/figs/network-035a-", paste(names(data)[blocks[, 
    j]], collapse = "-"), "-", gsub("\\.", "_", as.character(cutoff)))
nfn

## [1] "../reports/figs/network-035a-proteomics-metabolomics-0_7"

png(paste0(nfn, ".png"), res = 600, width = 4000, height = 4000)
nw <- my.network(final.diablo.model, blocks = blocks[, j], color.node = cols[blocks[, 
    j]], cutoff = cutoff, shape.node = "rectangle", save = "png", name.save = nfn)

Cutoff = 0.7

../reports/figs/network-035a-proteomics-metabolomics-0_7

{r echo=FALSE}gsub("\\.","_",as.character(cutoff)############################################### ## ## ## (c) Andrej Blejec (andrej.blejec@nib.si) ## ## ## ###############################################

3.6.2 proteomics and hormonomics

nfn <- paste0("../reports/figs/network-035a-", paste(names(data)[blocks[, 
    j]], collapse = "-"), "-", gsub("\\.", "_", as.character(cutoff)))
nfn

## [1] "../reports/figs/network-035a-proteomics-hormonomics-0_7"

png(paste0(nfn, ".png"), res = 600, width = 4000, height = 4000)
nw <- my.network(final.diablo.model, blocks = blocks[, j], color.node = cols[blocks[, 
    j]], cutoff = cutoff, shape.node = "rectangle", save = "png", name.save = nfn)

Cutoff = 0.7

../reports/figs/network-035a-proteomics-hormonomics-0_7

{r echo=FALSE}gsub("\\.","_",as.character(cutoff)############################################### ## ## ## (c) Andrej Blejec (andrej.blejec@nib.si) ## ## ## ###############################################

3.6.3 proteomics and qPCR

nfn <- paste0("../reports/figs/network-035a-", paste(names(data)[blocks[, 
    j]], collapse = "-"), "-", gsub("\\.", "_", as.character(cutoff)))
nfn

## [1] "../reports/figs/network-035a-proteomics-qPCR-0_7"

png(paste0(nfn, ".png"), res = 600, width = 4000, height = 4000)
nw <- my.network(final.diablo.model, blocks = blocks[, j], color.node = cols[blocks[, 
    j]], cutoff = cutoff, shape.node = "rectangle", save = "png", name.save = nfn)

Cutoff = 0.7

../reports/figs/network-035a-proteomics-qPCR-0_7

{r echo=FALSE}gsub("\\.","_",as.character(cutoff)############################################### ## ## ## (c) Andrej Blejec (andrej.blejec@nib.si) ## ## ## ###############################################

3.6.4 proteomics and phenomics

nfn <- paste0("../reports/figs/network-035a-", paste(names(data)[blocks[, 
    j]], collapse = "-"), "-", gsub("\\.", "_", as.character(cutoff)))
nfn

## [1] "../reports/figs/network-035a-proteomics-phenomics-0_7"

png(paste0(nfn, ".png"), res = 600, width = 4000, height = 4000)
nw <- my.network(final.diablo.model, blocks = blocks[, j], color.node = cols[blocks[, 
    j]], cutoff = cutoff, shape.node = "rectangle", save = "png", name.save = nfn)

Cutoff = 0.7

../reports/figs/network-035a-proteomics-phenomics-0_7

{r echo=FALSE}gsub("\\.","_",as.character(cutoff)############################################### ## ## ## (c) Andrej Blejec (andrej.blejec@nib.si) ## ## ## ###############################################

3.6.5 metabolomics and hormonomics

nfn <- paste0("../reports/figs/network-035a-", paste(names(data)[blocks[, 
    j]], collapse = "-"), "-", gsub("\\.", "_", as.character(cutoff)))
nfn

## [1] "../reports/figs/network-035a-metabolomics-hormonomics-0_7"

png(paste0(nfn, ".png"), res = 600, width = 4000, height = 4000)
nw <- my.network(final.diablo.model, blocks = blocks[, j], color.node = cols[blocks[, 
    j]], cutoff = cutoff, shape.node = "rectangle", save = "png", name.save = nfn)

Cutoff = 0.7

../reports/figs/network-035a-metabolomics-hormonomics-0_7

{r echo=FALSE}gsub("\\.","_",as.character(cutoff)############################################### ## ## ## (c) Andrej Blejec (andrej.blejec@nib.si) ## ## ## ###############################################

3.6.6 metabolomics and qPCR

nfn <- paste0("../reports/figs/network-035a-", paste(names(data)[blocks[, 
    j]], collapse = "-"), "-", gsub("\\.", "_", as.character(cutoff)))
nfn

## [1] "../reports/figs/network-035a-metabolomics-qPCR-0_7"

png(paste0(nfn, ".png"), res = 600, width = 4000, height = 4000)
nw <- my.network(final.diablo.model, blocks = blocks[, j], color.node = cols[blocks[, 
    j]], cutoff = cutoff, shape.node = "rectangle", save = "png", name.save = nfn)

Cutoff = 0.7

../reports/figs/network-035a-metabolomics-qPCR-0_7

{r echo=FALSE}gsub("\\.","_",as.character(cutoff)############################################### ## ## ## (c) Andrej Blejec (andrej.blejec@nib.si) ## ## ## ###############################################

3.6.7 metabolomics and phenomics

nfn <- paste0("../reports/figs/network-035a-", paste(names(data)[blocks[, 
    j]], collapse = "-"), "-", gsub("\\.", "_", as.character(cutoff)))
nfn

## [1] "../reports/figs/network-035a-metabolomics-phenomics-0_7"

png(paste0(nfn, ".png"), res = 600, width = 4000, height = 4000)
nw <- my.network(final.diablo.model, blocks = blocks[, j], color.node = cols[blocks[, 
    j]], cutoff = cutoff, shape.node = "rectangle", save = "png", name.save = nfn)

Cutoff = 0.7

../reports/figs/network-035a-metabolomics-phenomics-0_7

{r echo=FALSE}gsub("\\.","_",as.character(cutoff)############################################### ## ## ## (c) Andrej Blejec (andrej.blejec@nib.si) ## ## ## ###############################################

3.6.8 hormonomics and qPCR

nfn <- paste0("../reports/figs/network-035a-", paste(names(data)[blocks[, 
    j]], collapse = "-"), "-", gsub("\\.", "_", as.character(cutoff)))
nfn

## [1] "../reports/figs/network-035a-hormonomics-qPCR-0_7"

png(paste0(nfn, ".png"), res = 600, width = 4000, height = 4000)
nw <- my.network(final.diablo.model, blocks = blocks[, j], color.node = cols[blocks[, 
    j]], cutoff = cutoff, shape.node = "rectangle", save = "png", name.save = nfn)

Cutoff = 0.7

../reports/figs/network-035a-hormonomics-qPCR-0_7

{r echo=FALSE}gsub("\\.","_",as.character(cutoff)############################################### ## ## ## (c) Andrej Blejec (andrej.blejec@nib.si) ## ## ## ###############################################

3.6.9 hormonomics and phenomics

nfn <- paste0("../reports/figs/network-035a-", paste(names(data)[blocks[, 
    j]], collapse = "-"), "-", gsub("\\.", "_", as.character(cutoff)))
nfn

## [1] "../reports/figs/network-035a-hormonomics-phenomics-0_7"

png(paste0(nfn, ".png"), res = 600, width = 4000, height = 4000)
nw <- my.network(final.diablo.model, blocks = blocks[, j], color.node = cols[blocks[, 
    j]], cutoff = cutoff, shape.node = "rectangle", save = "png", name.save = nfn)

Cutoff = 0.7

../reports/figs/network-035a-hormonomics-phenomics-0_7

{r echo=FALSE}gsub("\\.","_",as.character(cutoff)############################################### ## ## ## (c) Andrej Blejec (andrej.blejec@nib.si) ## ## ## ###############################################

3.6.10 qPCR and phenomics

nfn <- paste0("../reports/figs/network-035a-", paste(names(data)[blocks[, 
    j]], collapse = "-"), "-", gsub("\\.", "_", as.character(cutoff)))
nfn

## [1] "../reports/figs/network-035a-qPCR-phenomics-0_7"

png(paste0(nfn, ".png"), res = 600, width = 4000, height = 4000)
nw <- my.network(final.diablo.model, blocks = blocks[, j], color.node = cols[blocks[, 
    j]], cutoff = cutoff, shape.node = "rectangle", save = "png", name.save = nfn)

Cutoff = 0.7

../reports/figs/network-035a-qPCR-phenomics-0_7

3.7 Multipartite network

cutoff <- 0
x <- final.diablo.model
layout.fun <- NULL
label <- paste(.treat, collapse = ", ")
out35b <- ""
out35b <- paste(out35b, knit_child("035b-multipartite-network.Rmd", quiet = TRUE))
cat(out35b)

3.7.1 Cutoff = 0

ndata <- length(data)
lbl <- gsub(", ", "-", label)
nfn <- paste("../reports/figs/network-035b", lbl, gsub("\\.", "_", as.character(cutoff)), 
    sep = "-")

set.seed(1234)
nw <- my.network(x, blocks = 1:ndata, color.node = cols, cutoff = cutoff, 
    shape.node = "rectangle", layout = layout.fun, save = "png", name.save = nfn)

../reports/figs/network-035b-H-HD-0

# Save network layout in ly, used by my.layout function.
if (exists(deparse(substitute(nw)))) ly <- nw$layout else ly <- NULL

cutoff <- 0.7
x <- final.diablo.model
label <- paste(.treat, collapse = ", ")
out35b <- ""
out35b <- paste(out35b, knit_child("035b-multipartite-network.Rmd", quiet = TRUE))
cat(out35b)

3.7.2 Cutoff = 0.7

ndata <- length(data)
lbl <- gsub(", ", "-", label)
nfn <- paste("../reports/figs/network-035b", lbl, gsub("\\.", "_", as.character(cutoff)), 
    sep = "-")

set.seed(1234)
nw <- my.network(x, blocks = 1:ndata, color.node = cols, cutoff = cutoff, 
    shape.node = "rectangle", layout = layout.fun, save = "png", name.save = nfn)

../reports/figs/network-035b-H-HD-0_7

3.8 Visualize loadings

The function “plotLoadings” visualises the loading weights of each selected variable on each component and each data set. The colour indicates the class in which the variable has the maximum level of expression “contrib = ‘max’” using the median “method = ‘median’”. Figure below depicts the loading values for each dimension.

for (i in 1:ncomp) plotLoadings(final.diablo.model, comp = i, contrib = "max", 
    method = "median")

3.9 Clustered image map ( Heatmap )

The cimDIABLO() function is a clustered image map specifically implemented to represent the multi-omics molecular signature expression for each sample. From figure below the areas of homogeneous expression levels for a set of samples across a set of features can be determined. For instance, the H14 samples were the only group to show extremely high levels of expression for a specific set of genes and metabolites. This indicates these features are fairly discriminating for this subtype.

cimfn <- "../reports/figs/cim.png"

png(cimfn, res = 600, width = 4000, height = 4000)
cim <- cimDiablo(final.diablo.model, size.legend = 0.7, comp = 1:3)
dev.off()

## pdf 
##   2

../reports/figs/cim.png

3.10 AUC and ROC plots

An AUC plot per block can also be obtained using the function auroc(). The interpretation of this output may not be particularly insightful in relation to the performance evaluation of our methods, but can complement the statistical analysis.

par(mfrow = c(2, 2))
for (i in 1:length(data)) auc.splsda = auroc(final.diablo.model, roc.block = names(data[i]), 
    roc.comp = 1, print = FALSE)

Save finil DIABLO model for future use in networks.

res <- final.diablo.model

Estimate classification error rate. The error rate should drop by more components used.

# run component number tuning with repeated CV
system.time(perf.diablo <- perf(res, validation = "Mfold", folds = 3, nrepeat = 10))

##    user  system elapsed 
##   14.57    0.33   15.03

plot(perf.diablo)  # plot output of tuning

3.11 Names of kept variables

# the features selected to form components
for (comp in 1:ncomp) {
    cat("\nComponent ", comp, ":\n")
    for (i in 1:length(data)) {
        cat(names(data)[i], "\n")
        print(selectVar(res, comp = comp)[[i]]$name)
    }
}

## 
## Component  1 :
## proteomics 
## [1] "PBdnRY1_2239"    "PBdnRY1_5389"    "PBdnRY1_7025"   
## [4] "CLCdnDe6_4953"   "CLCdnPW49_8560"  "VdnPW3_37011"   
## [7] "SdnPW1_46"       "CLCdnPW10_11437"
## metabolomics 
## [1] "Glukose"  "Fructose" "Ile"      "Starch"   "Val"      "Gly"     
## [7] "His"      "Tyr"     
## hormonomics 
## [1] "SA"         "ABA"        "X12.OH.JA"  "DPA"        "JA.Ile"    
## [6] "PA"         "cisOPDA"    "X9.10.dhJA"
## qPCR 
## [1] "RD29B"   "X13.LOX" "SP6A"    "M0ZJG3"  "RbohA"   "HSP70"  
## [7] "CAT1"    "PR1b"   
## phenomics 
## [1] "Fv.Fm"             "DW_SH_Total"       "deltaT"           
## [4] "qL"                "FW_SH_Total"       "TOP.AREA"         
## [7] "COMPACTNESS"       "WATER.CONSUMPTION"
## 
## Component  2 :
## proteomics 
## [1] "VdnDe2_53978"       "CLCdnPW48_15164"    "CLCdnPW12_1472"    
## [4] "Sotub11g025210.1.1" "CLCdnRY10_6741"     "CLCdnPW31_956"     
## [7] "VdnPW4_15664"       "CLCdnPW22_15452"   
## metabolomics 
## [1] "Ala" "Gln" "Met" "Asn" "Thr" "Leu" "Gly" "Arg"
## hormonomics 
## [1] "oxIAA"      "JA"         "X9.10.dhJA" "IAA.Asp"    "ABA"       
## [6] "IAA"        "X12.OH.JA"  "DPA"       
## qPCR 
## [1] "ACO2"  "SnRK2" "ERF1"  "PR1b"  "RbohA" "HSP70" "CAT1"  "P5CS" 
## phenomics 
## [1] "TOP.AREA"          "WATER.CONSUMPTION" "COMPACTNESS"      
## [4] "Fv.Fm"             "FW_SH_Total"       "deltaT"           
## [7] "DW_SH_Total"       "qL"               
## 
## Component  3 :
## proteomics 
## [1] "CLCdnDe7_6144"   "SdnRY1_10712"    "VdnDe1_39992"   
## [4] "SdnPW1_46"       "CLCdnPW49_8560"  "CLCdnPW31_956"  
## [7] "VdnPW4_18801"    "CLCdnDe13_63179"
## metabolomics 
## [1] "Asp"     "Leu"     "Lys"     "Ser"     "Pro"     "His"    
## [7] "Sucrose" "Arg"    
## hormonomics 
## [1] "X12.OH.JA"  "X9.10.dhJA" "cisOPDA"    "SA"         "PA"        
## [6] "JA"         "JA.Ile"     "DPA"       
## qPCR 
## [1] "RbohA"   "PR1b"    "RD29B"   "ERF1"    "P5CS"    "SP6A"   
## [7] "ACO2"    "X13.LOX"
## phenomics 
## [1] "COMPACTNESS"       "qL"                "deltaT"           
## [4] "WATER.CONSUMPTION" "DW_SH_Total"       "FW_SH_Total"      
## [7] "Fv.Fm"             "TOP.AREA"         
## 
## Component  4 :
## proteomics 
## [1] "CLCdnPW10_11437" "VdnDe4_115100"   "VdnPW4_18801"   
## [4] "VdnPW5_1541"     "CLCdnPW22_15452" "VdnDe2_86784"   
## [7] "CLCdnDe10_24187" "CLCdnPW12_1472" 
## metabolomics 
## [1] "Arg"     "Phe"     "His"     "Asn"     "Asp"     "Met"    
## [7] "Sucrose" "Tyr"    
## hormonomics 
## [1] "IAA"        "IAA.Asp"    "JA.Ile"     "PA"         "X9.10.dhJA"
## [6] "JA"         "X12.OH.JA"  "ABA"       
## qPCR 
## [1] "P5CS"    "ERF1"    "M0ZJG3"  "PR1b"    "SP6A"    "X13.LOX"
## [7] "RD29B"   "RbohA"  
## phenomics 
## [1] "WATER.CONSUMPTION" "Fv.Fm"             "FW_SH_Total"      
## [4] "DW_SH_Total"       "qL"                "deltaT"           
## [7] "COMPACTNESS"       "TOP.AREA"         
## 
## Component  5 :
## proteomics 
## [1] "VdnDe2_86784"           "TDe_comp110944_c2_seq4"
## [3] "VdnPW5_1541"            "Sotub11g025210.1.1"    
## [5] "CLCdnDe2_311"           "PBdnRY1_2239"          
## [7] "VdnDe1_162036"          "CLCdnPW22_15452"       
## metabolomics 
## [1] "Glu" "Asn" "Leu" "Ser" "His" "Pro" "Tyr" "Thr"
## hormonomics 
## [1] "PA"        "JA.Ile"    "cisOPDA"   "IAA"       "X12.OH.JA"
## [6] "SA"        "JA"        "ABA"      
## qPCR 
## [1] "PR1b"    "SP6A"    "X13.LOX" "RD29B"   "CAT1"    "ERF1"   
## [7] "SWEET"   "M0ZJG3" 
## phenomics 
## [1] "qL"                "DW_SH_Total"       "COMPACTNESS"      
## [4] "WATER.CONSUMPTION" "deltaT"            "FW_SH_Total"      
## [7] "TOP.AREA"          "Fv.Fm"

One would like to reduce the number of nodes, especially for proteomics data. One option is to reduce datasets in a way to keep only the variables in the selectVars in original data in . We will keep variables from the first two components.

keptVars <- unique(c(selectVar(res, comp = 1)[[1]]$name, selectVar(res, 
    comp = 2)[[1]]$name))
which(keptVars %in% selectVar(res, comp = 1)[[1]]$name)

## [1] 1 2 3 4 5 6 7 8

which(keptVars %in% selectVar(res, comp = 2)[[1]]$name)

## [1]  9 10 11 12 13 14 15 16

3.12 Networks circos plots

size.variables <- 1
sim <- circosPlot(final.diablo.model, cutoff = 0.5, line = TRUE, color.blocks = cols, 
    color.cor = c(3, 2), size.labels = 1, size.variables = size.variables, 
    xpd = TRUE)

circosPlot(final.diablo.model, cutoff = 0.75, line = TRUE, color.blocks = cols, 
    color.cor = c(3, 2), size.labels = 1, size.variables = size.variables, 
    xpd = TRUE)

circosPlot(final.diablo.model, cutoff = 0.9, line = TRUE, color.blocks = cols, 
    color.cor = c(3, 2), size.labels = 1, size.variables = size.variables, 
    xpd = TRUE)

4 Networks

We will prepare partial models for each treatment and treatment combination.

netlist <- list()
layout.fun <- NULL
for (i in 1:length(.treat)) {
    thisTreat <- .treat[i]
    out037 <- ""
    out037 <- paste(out037, knit_child("037-network.Rmd", quiet = TRUE))
    netlist[[thisTreat]] <- RETURN  # complete network
    
    cat(out037)
}

4.1 Network for H.

filter <- pdata$Treatment %in% thisTreat
XX1 <- lapply(CCDATA, function(x) if (is.null(dim(x))) x[filter] else x[filter, 
    ])
table(XX1$status)

## 
## H01 H07 H08 H14 
##   4   4   4   4

res1 <- block.splsda(X = XX1[-1], Y = as.factor(XX1[[1]]), ncomp = ncomp, 
    keepX = list.keepX, design = design)

## Design matrix has changed to include Y; each block will be
##             linked to Y.

cutoff <- 0
x <- res1
layout.fun <- NULL
label <- thisTreat
out23b <- ""
out23b <- paste(out23b, knit_child("035b-multipartite-network.Rmd", quiet = TRUE))
RETURN <- nw  # complete network
cat(out23b)

4.1.1 Cutoff = 0

ndata <- length(data)
lbl <- gsub(", ", "-", label)
nfn <- paste("../reports/figs/network-035b", lbl, gsub("\\.", "_", as.character(cutoff)), 
    sep = "-")

set.seed(1234)
nw <- my.network(x, blocks = 1:ndata, color.node = cols, cutoff = cutoff, 
    shape.node = "rectangle", layout = layout.fun, save = "png", name.save = nfn)

../reports/figs/network-035b-H-0

Save network layout for further plots, used by layout function my.layout.

ly <- nw$layout

cutoff <- 0.7
x <- res1
layout.fun <- my.layout
label <- thisTreat
out23b <- ""
out23b <- paste(out23b, knit_child("035b-multipartite-network.Rmd", quiet = TRUE))
cat(out23b)

4.1.2 Cutoff = 0.7

ndata <- length(data)
lbl <- gsub(", ", "-", label)
nfn <- paste("../reports/figs/network-035b", lbl, gsub("\\.", "_", as.character(cutoff)), 
    sep = "-")

set.seed(1234)
nw <- my.network(x, blocks = 1:ndata, color.node = cols, cutoff = cutoff, 
    shape.node = "rectangle", layout = layout.fun, save = "png", name.save = nfn)

../reports/figs/network-035b-H-0_7

4.2 Network for HD.

filter <- pdata$Treatment %in% thisTreat
XX1 <- lapply(CCDATA, function(x) if (is.null(dim(x))) x[filter] else x[filter, 
    ])
table(XX1$status)

## 
## HD08 HD14 
##    4    4

res1 <- block.splsda(X = XX1[-1], Y = as.factor(XX1[[1]]), ncomp = ncomp, 
    keepX = list.keepX, design = design)

## Design matrix has changed to include Y; each block will be
##             linked to Y.

cutoff <- 0
x <- res1
layout.fun <- NULL
label <- thisTreat
out23b <- ""
out23b <- paste(out23b, knit_child("035b-multipartite-network.Rmd", quiet = TRUE))
RETURN <- nw  # complete network
cat(out23b)

4.2.1 Cutoff = 0

ndata <- length(data)
lbl <- gsub(", ", "-", label)
nfn <- paste("../reports/figs/network-035b", lbl, gsub("\\.", "_", as.character(cutoff)), 
    sep = "-")

set.seed(1234)
nw <- my.network(x, blocks = 1:ndata, color.node = cols, cutoff = cutoff, 
    shape.node = "rectangle", layout = layout.fun, save = "png", name.save = nfn)

../reports/figs/network-035b-HD-0

Save network layout for further plots, used by layout function my.layout.

ly <- nw$layout

cutoff <- 0.7
x <- res1
layout.fun <- my.layout
label <- thisTreat
out23b <- ""
out23b <- paste(out23b, knit_child("035b-multipartite-network.Rmd", quiet = TRUE))
cat(out23b)

4.2.2 Cutoff = 0.7

ndata <- length(data)
lbl <- gsub(", ", "-", label)
nfn <- paste("../reports/figs/network-035b", lbl, gsub("\\.", "_", as.character(cutoff)), 
    sep = "-")

set.seed(1234)
nw <- my.network(x, blocks = 1:ndata, color.node = cols, cutoff = cutoff, 
    shape.node = "rectangle", layout = layout.fun, save = "png", name.save = nfn)

../reports/figs/network-035b-HD-0_7

names(netlist)

## [1] "H"  "HD"

4.3 Network for H and HD.

# Complete network, cutoff = 0, all treatments
datasets <- names(CCDATA[-1])
ndatasets <- length(datasets)
# 
nname <- paste("../ouptut/figs/network-", paste(useTreatment, collapse = ""), 
    sep = "")
N12 <- my.network(res, cutoff = 0, blocks = 1:ndatasets, shape.node = c("rectangle"), 
    save = "png", name.save = "network-CH")
# 
e <- extractEdges2(N12)
colnames(e)[ncol(e)] <- paste(.treat, collapse = ".")
head(e)

##                                                        edge
## pr.VdnDe2_86784_me.Glukose       pr.VdnDe2_86784_me.Glukose
## pr.PBdnRY1_7025_me.Glukose       pr.PBdnRY1_7025_me.Glukose
## pr.VdnPW5_1541_me.Glukose         pr.VdnPW5_1541_me.Glukose
## pr.CLCdnPW10_11437_me.Glukose pr.CLCdnPW10_11437_me.Glukose
## pr.CLCdnDe6_4953_me.Glukose     pr.CLCdnDe6_4953_me.Glukose
## pr.CLCdnPW12_1472_me.Glukose   pr.CLCdnPW12_1472_me.Glukose
##                                   group1            from       group2
## pr.VdnDe2_86784_me.Glukose    proteomics    VdnDe2_86784 metabolomics
## pr.PBdnRY1_7025_me.Glukose    proteomics    PBdnRY1_7025 metabolomics
## pr.VdnPW5_1541_me.Glukose     proteomics     VdnPW5_1541 metabolomics
## pr.CLCdnPW10_11437_me.Glukose proteomics CLCdnPW10_11437 metabolomics
## pr.CLCdnDe6_4953_me.Glukose   proteomics   CLCdnDe6_4953 metabolomics
## pr.CLCdnPW12_1472_me.Glukose  proteomics  CLCdnPW12_1472 metabolomics
##                                    to       H.HD
## pr.VdnDe2_86784_me.Glukose    Glukose -0.3679799
## pr.PBdnRY1_7025_me.Glukose    Glukose -0.7343140
## pr.VdnPW5_1541_me.Glukose     Glukose -0.3003502
## pr.CLCdnPW10_11437_me.Glukose Glukose -0.5412361
## pr.CLCdnDe6_4953_me.Glukose   Glukose -0.8443068
## pr.CLCdnPW12_1472_me.Glukose  Glukose  0.5979807

tail(e)

##                                                          edge group1
## qP.P5CS_ph.WATER.CONSUMPTION     qP.P5CS_ph.WATER.CONSUMPTION   qPCR
## qP.ERF1_ph.WATER.CONSUMPTION     qP.ERF1_ph.WATER.CONSUMPTION   qPCR
## qP.CAT1_ph.WATER.CONSUMPTION     qP.CAT1_ph.WATER.CONSUMPTION   qPCR
## qP.SWEET_ph.WATER.CONSUMPTION   qP.SWEET_ph.WATER.CONSUMPTION   qPCR
## qP.SP6A_ph.WATER.CONSUMPTION     qP.SP6A_ph.WATER.CONSUMPTION   qPCR
## qP.M0ZJG3_ph.WATER.CONSUMPTION qP.M0ZJG3_ph.WATER.CONSUMPTION   qPCR
##                                  from    group2                to
## qP.P5CS_ph.WATER.CONSUMPTION     P5CS phenomics WATER.CONSUMPTION
## qP.ERF1_ph.WATER.CONSUMPTION     ERF1 phenomics WATER.CONSUMPTION
## qP.CAT1_ph.WATER.CONSUMPTION     CAT1 phenomics WATER.CONSUMPTION
## qP.SWEET_ph.WATER.CONSUMPTION   SWEET phenomics WATER.CONSUMPTION
## qP.SP6A_ph.WATER.CONSUMPTION     SP6A phenomics WATER.CONSUMPTION
## qP.M0ZJG3_ph.WATER.CONSUMPTION M0ZJG3 phenomics WATER.CONSUMPTION
##                                      H.HD
## qP.P5CS_ph.WATER.CONSUMPTION   0.18115089
## qP.ERF1_ph.WATER.CONSUMPTION   0.31092616
## qP.CAT1_ph.WATER.CONSUMPTION   0.17666551
## qP.SWEET_ph.WATER.CONSUMPTION  0.16032835
## qP.SP6A_ph.WATER.CONSUMPTION   0.02235330
## qP.M0ZJG3_ph.WATER.CONSUMPTION 0.09472415

dim(e)

## [1] 2622    6

Save networks file for combined and single treatments, stored in netlist. Combine edges from all networks as columns in a data.frame.

for (i in 1:length(netlist)) {
    e1 <- extractEdges2(netlist[[i]])
    colnames(e1)[ncol(e1)] <- names(netlist)[i]
    e <- merge(e, e1, sort = FALSE, all = TRUE)
}

str(e)

## 'data.frame':    2908 obs. of  8 variables:
##  $ edge  : chr  "pr.PBdnRY1_7025_me.Glukose" "pr.VdnPW5_1541_me.Glukose" "pr.CLCdnPW10_11437_me.Glukose" "pr.CLCdnDe6_4953_me.Glukose" ...
##  $ group1: chr  "proteomics" "proteomics" "proteomics" "proteomics" ...
##  $ from  : chr  "PBdnRY1_7025" "VdnPW5_1541" "CLCdnPW10_11437" "CLCdnDe6_4953" ...
##  $ group2: chr  "metabolomics" "metabolomics" "metabolomics" "metabolomics" ...
##  $ to    : chr  "Glukose" "Glukose" "Glukose" "Glukose" ...
##  $ H.HD  : num  -0.734 -0.3 -0.541 -0.844 0.598 ...
##  $ H     : num  -0.514 -0.174 -0.671 -0.651 0.383 ...
##  $ HD    : num  -0.748 -0.386 -0.444 -0.867 0.743 ...

Compose file name and necessary information for network export file

netfn <- paste0("../output/network-", paste(.treat, collapse = "_"), "-", 
    paste(datasets, collapse = "_"), ".txt")
label0 <- paste(paste(.treat, collapse = ", "), "|", paste(datasets, collapse = ", "), 
    "; cutoff =", 0)
title <- label0
sets <- 1:length(DATA)
suffix <- paste0(substr(names(DATA), 1, 2)[sets[-1]], collapse = "-")

4.4 Export edges table

# e <- data.frame(x=1:10,y=1:10) my.write.table(e,
# file='network.txt',meta=FALSE)
write.table(e, file = netfn, na = "0", sep = "\t", col.names = NA, row.names = TRUE)

Network edges are exported into file

netfn

## [1] "../output/network-H_HD-proteomics_metabolomics_hormonomics_qPCR_phenomics.txt"

Table with edges for networks based on combined treatments (C, H) and single treatments (C) and (H) is exported as a text file. This table can be used for inspection and filtering out edges based on selected cutoff. Missing edges are labeled as weight 0. This enables numeric filtration in other visualization or analysis files (e.g. Excel).

5 SessionInfo

Windows 10 x64 (build 19045)
R version 4.0.2 (2020-06-22) Platform: x86_64-w64-mingw32/x64 (64-bit) Running under: Windows 10 x64 (build 19045)

Matrix products: default

locale: [1] LC_COLLATE=Slovenian_Slovenia.1250 [2] LC_CTYPE=Slovenian_Slovenia.1250
[3] LC_MONETARY=Slovenian_Slovenia.1250 [4] LC_NUMERIC=C
[5] LC_TIME=Slovenian_Slovenia.1250
system code page: 1252

attached base packages: [1] grid stats graphics utils datasets grDevices [7] methods base

other attached packages: [1] pheatmap_1.0.12 ComplexHeatmap_2.6.2 igraph_1.2.6
[4] mixOmics_6.14.0 ggplot2_3.3.5 lattice_0.22-5
[7] MASS_7.3-60.0.1 knitr_1.43

loaded via a namespace (and not attached): [1] ggrepel_0.9.0 Rcpp_1.0.7 circlize_0.4.15
[4] tidyr_1.1.2 corpcor_1.6.9 png_0.1-7
[7] digest_0.6.35 RSpectra_0.16-0 R6_2.5.1
[10] plyr_1.8.6 stats4_4.0.2 ellipse_0.4.2
[13] evaluate_0.21 highr_0.8 pillar_1.4.7
[16] GlobalOptions_0.1.2 rlang_1.1.1 rstudioapi_0.13
[19] jquerylib_0.1.4 S4Vectors_0.28.1 GetoptLong_1.0.5
[22] Matrix_1.6-5 rmarkdown_2.21 labeling_0.4.2
[25] rARPACK_0.11-0 stringr_1.4.0 munsell_0.5.0
[28] compiler_4.0.2 xfun_0.39 pkgconfig_2.0.3
[31] BiocGenerics_0.36.0 shape_1.4.6 htmltools_0.5.2
[34] tidyselect_1.1.0 tibble_3.0.4 gridExtra_2.3
[37] IRanges_2.24.1 matrixStats_1.2.0 crayon_1.3.4
[40] dplyr_1.0.2 withr_3.0.0 jsonlite_1.8.8
[43] gtable_0.3.0 lifecycle_0.2.0 magrittr_2.0.1
[46] formatR_1.7 scales_1.1.1 stringi_1.5.3
[49] farver_2.0.3 reshape2_1.4.4 bslib_0.3.1
[52] ellipsis_0.3.1 generics_0.1.0 vctrs_0.4.1
[55] rjson_0.2.20 RColorBrewer_1.1-2 tools_4.0.2
[58] Cairo_1.5-15 glue_1.4.2 purrr_0.3.4
[61] parallel_4.0.2 fastmap_1.1.1 yaml_2.2.1
[64] clue_0.3-60 colorspace_1.4-1 cluster_2.1.0
[67] sass_0.4.0